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FLOW  IN  COMBUSTION  SYSTEMS 


Abstract 


Combustion  systems  involve  the  complex  interaction  between  several 
fundamental  phenomena^hence  the  fundamental  aspects  of  each  phenomena  cannot 
be  studied  in  isolation  if  we  wish  to  understiandnEfie^whole  processT"  Vln 
this  investigation,  the  basic  science  underlying  the  interactions  between  the 
two-phase  flow,  fluid  dynamics  and  chemical  kinetics  have  been  investigated. 
The  studies  have  required  the  development  of  new  diagnostic  systems  and 
significant  progress  has  been  made  in  the  following  areas :- 


The  development  of  a  technique  for  making  accurate  dropsize  measurements 
in  dense  sprays. 


The  application  of  this  technique  to  an  F101  air  blast  atomizer. 

The  use  of  LDA  for  the  precise  characterization  of  swirl  from  the  F101 
swirler . 

The  development  of  shear  stress  mathematical  models  for  non-isotropic 
turbulence . 

The  application  of  this  model  to  the  FlOl  swirler. 

The  development  of  a  mercury  vapour  pulse  tracer  for  residence  time 
distribution  measurement  in  combustors. 


The  development  of  a  mathematical  modelling  technique  whereby  the 
residence  time  distribution  can  be  computed. 

Closing  of  the  gap  between  stirred  reactor  models  and  finite  difference 
models  of  combustion  systems. 

— Proposal  of  a  new  fundamental  approach  to  the  problem  of  simultaneous 
mixing  and  reaction  using  a  quantitiative  coalescence/dispersion  eddy 
concept  which  has  the  potential  to  represent  all  the  high  order 
correlations  of  the  interaction.  Of, Ultluit'. 
STUDY  OF  THE  FlOl  ANNULAR  GAS  TURBINE  COMBUSTOR. 


The  object  of  this  study  is  to  apply  the  Sheffield  3-D  flow  model  to 
the  geometry  of  the  FlOl  combustion  chamber.  This  work,  together  with 
experimental  studies  carried  out  on  a  54°  sector  FlOl  chamber,  will  provide 
a  useful  testing  ground  for  future  model  development. 

To  this  aim,  a  calibration  rig  has  been  constructed  to  make  inlet 
condition  measurements  on  the  FlOl  combustor  spectacle  plate.  This  has  been 
used  to  examine  the  fuel  spray  droplet  size  distribution  and  to  study  the  FlOl 
counter  swirling  primary  air  inlet.  L.D.A.  measurements  in  the  swirling  flow 
have  been  compared  with  computer  predictions  and  results  are  encouraging. 

A  combustor  test  facility  has  been  installed  to  allow  cold  and  hot  flow 
measurements  to  be  made  both  at  atmospheric  and  elevated  pressures .  This 
facility  has  been  used  to  obtain  residence  time  distribution  measurements  using 
a  pulsed  mercury  tracer. 


The  work  is  described  in  detail  in  the  following  three  reports: 

Fuel  spray  droplet  sizing  on  the  F101  air  blast  atomizer. 

A  study  of  the  F101  combustor  swirlers. 

The  F101  combustor  test  facility  and  Mercury  pulse  tracer 
measurements . 

Other  relevant  papers  included  are 

P.G.  Felton,  A. A.  Hamidi  and  A.K.  Aigal  "Multiple  scattering  effects  on 
particle  size  by  laser  diffraction" . 

B.C.R.  Ewan,  F.  Boysan  and  J.  Swithenbank  "Closing  the  gap  between  finite 
difference  and  stirred  reactor  combustor  modelling". 


FUEL  SPRAY  DROPLET  SIZING  ON  THE  F101  AIR  BLAST  ATOMIZER. 


Introduction. 


Liquid  fuel  injected  into  gas  turbine  combustors  is  usually  atomized 
by  means  of  a  pressure  jet  or  air  blast  atomizer.  The  spray  produced  will 
consist  of  droplets  with  a  wide  range  of  diameters.  Achieving  the  correct 
droplet  size  distribution  over  all  operating  conditions  is  important  in 
combustor  design  since  it  determines  the  distribution  of  evaporated  fuel  in  the 
primary  zone.  Too  fine  a  droplet  will  not  penetrate  far  into  the  air  stream 
and  too  large  a  droplet  will  not  evaporate  rapidly  enough.  Earlier  work  at 
Sheffield  clearly  demonstrated  the  significance  of  these  effects. 


The  FLUENT  code  (1)  for  three  dimensional  two  phase  reacting  flows  is 
able  to  calculate  the  trajectories  of  individual  fuel  droplets  in  a  combustor 
and  to  simulate  their  effect  on  evaporated  fuel  distribution.  It  is  also  able 


to  predict  the  effect  of  turbulence  on  these  droplet  trajectories.  To  do  this 
the  initial  size  distribution  of  the  spray  is  divided  into  a  number  of  discrete 
intervals  each  representing  an  average  diameter.  The  equations  for  motion, 
heating  up  and  evaporation  are  then  solved  for  droplets  representing  each  size 
group. 

It  is  desired  to  model  the  F101  combustor  using  the  FLUENT  code  and  so 
the  fuel  droplet  size  distribution,  produced  by  its  air  blast  atomizer,  must  be 
determined.  This  report  described  how  this  is  achieved  using  a  Malvern 
particle  sizing  instrument. 


The  Malvern  Particle  Sizer,  Operating  Theory. 


The  Malvern  particle  sizer  provides  a  non-intrusive  optical  technique  for 
particle  sizing.  The  optical  configuration  is  shown  in  figure  5.  A  parallel 
beam  of  monochromatic  coherent  light,  produced  by  a  1  mW  He/Ne  laser  and  beam 
expander,  is  passed  through  the  field  of  droplets  or  particles  to  be  measured. 
When  a  spherical  particle  is  illuminated  in  this  manner  a  Fraunhofer  diffraction 
pattern  is  formed.  If  the  particle  is  arranged  to  be  within  the  focal  length  of 
a  Fourier  transform  lens  then  undiffracted  light  is  focused  to  a  point  on  the 
optical  axis  and  diffracted  light  forms  a  far  field  Fraunhofer  pattern  of  rings 
around  the  central  point.  The  diffraction  pattern  formed  by  a  single  particle 
or  a  field  of  monosized  particles  is  shown  in  figure  6.  The  pattern  will  be 
stationary,  independent  of  whether  the  particles  are  in  motion,  since  light  is 
diffracted  at  an  angle  dependent  only  on  particle  size  and  will  give  the  same 
radial  displacement  in  the  focal  plane  irrespective  of  particle  position  in  the 
illuminating  beam. 

For  a  monosized  field  of  particles,  of  size  a,  the  light  energy 
contained  within  a  circle,  radius  s,  in  the  focal  plane  is  given  by 


The  light  energy  within  any  ring,  in  the  focal  plane,  bounded  by  radii  Sj. 
and  s„  is  therefore  given  by 
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where 


E  is  the  energy  incident  on  the  particles  and  is  proportional  to  their  cross 
sectional  area  and  their  number . 


C  7T  a  N 


is  a  constant  dependent  on  laser  power  and  N  is  the  number  of  particles, 
given  by 
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Hence 
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If  a  range  of  particle  sizes  are  present  then  the  light  energy,  falling  onto 
any  ring  in  the  focal  plane,  is  the  sum  of  the  contributions  from  all 
particles. 
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The  particle  size  distribution  is  classified  into  a  number  of  size 
ranges  and  incident  light  energy  is  measured  by  means  of  concentric  photo- 
detector  rings  placed  in  the  focal  plane.  Measurements  are  made  at  radii 


where  s  is  related  to  a^,  the  mean  size  of  the  ith  size  range,  by 
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Equation  (8)  is  the  condition  for  the  first  maximum  in  the  diffracted  energy 
distribution. 


The  total  energy  distribution  is  the  sum  of  the  products  of  the  energy 
distributions  for  each  size  range  and  the  weight  fraction  of  particles  in  that 
range.  In  matrix  form 
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Where  L(I)  is  the  light  energy  falling  on  ring  I;  W(J)  is  the  weight  fraction 
in  size  range  J  and  T(I,J)  contains  the  coefficients  which  define  the  light 
energy  distribution  for  each  particle  size  range. 

A  weight  distribution,  W(J) ,  is  fitted  to  the  measured  distribution  of 
light  energy  by  an  iterative  procedure  such  that  the  error 

l  (L (I)  -  W(J)T(I ,J) ) 2  (10) 


is  minimized. 

This  analysis  is  carried  out  by  a  dedicated  CBM  microcomputer.  The 
system  allows  any  of  a  number  of  different  size  distribution  functions  to  be 
fitted  to  the  measured  energy  distribution.  These  include  Rosin  Rammler 
and  log  normal  distribution  as  well  as  a  modal  independent  option. 

For  fuel  spray  studies  the  Rosin  Rammler  distribution  is  often  employed 
as  this  has  been  found  to  fit  spray  size  distributions  particularly  well. 

The  distribution  is  defined  by 

Y  =  10  s  ( 1-Exp (-x/x)N)  (11) 

Y  =  %  of  particles  with  diameter  less  than  X. 

X  =  Rosin  Rammler  mean  diameter. 

N  =  Rosin  Rammler  exponent  which  measures  the  spread  of  the 

distribution  about  the  mean  X. 


The  difference  between  the  light  intensity  measured  at  the  centre  of 
the  photodetector,  before  and  after  particles  are  placed  in  the  beam  gives 
the  fraction  of  light  obscured  (or  the  obscuration) .  This  is  related  to  the 
total  projected  area  of  the  particles  by  the  Beer  Lambert  law 
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and  allows  determination  of  particle  concentration. 
Limitations  of  the  Malvern  Sizing  Instrument. 


L 


The  instrument  is  capable  of  measuring  particle  size  down  to  one  micron. 
Smaller  particles  cannot  be  measured  since  the  light  scatterers  must  be  larger 
than  the  wavelength  of  light  (0.6238  um  for  a  He/Ne  laser). 
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All  the  particles  must  be  within  the  focal  length  of  the  Fourier  transform 
lens.  Fine  particles  which  are  too  far  from  the  lens  will  scatter  light  which 
is  not  focused  onto  the  photodetectors  and  results  will  be  biased  towards  large 
particles. 


Reliable  results  can  only  be  obtained  if  the  light  obscuration  is 
in  the  range  5  to  50%.  Low  obscuration  means  that  a  representative  sample 
of  particles  is  not  analysed  while  at  obscurations  over  50%  multiple 
scattering  of  light  becomes  significant  and  biases  results  towards  small 
particles.  Hamidi  and  Felton  (2)  have  recently  developed  a  model  for 
the  latter  effect  which  allows  corrections  for  multiple  scattering  to  be 
made. 

The  F101  Fuel  Injector. 

Figure  1(a)  shows  the  F101  CED  fuel  injector  nozzle.  Fuel  enters 
through  four  ports  in  the  nozzle  and  annular  openings  around  these  surround 
the  fuel  with  high  velocity  air.  The  nozzles  are  mounted  in  the  spectacle 
plate  shown  in  figure  1(b) .  Fuel  from  each  nozzle  is  carried  into  the  com¬ 
bustor  primary  zone  by  air  from  two  concentric  counter  swirling  inlets  in 
the  spectacle  plate. 

Sizing  Experiment. 

The  apparatus  used  for  this  is  illustrated  in  figure  2.  Measurements 
cannot  be  made  inside  the  combustor  since  geometry  limits  access.  The 
spectacle  plate  was  therefore  removed  and  mounted  on  a  milling  machine  table 
as  shown  in  the  figure.  This  provides  a  rigid  support  and  also  allows  accurate 
movement  of  the  plate  in  three  dimensions. 

Fuel  and  air  were  supplied  to  a  single  atomizer  nozzle.  Air  was  also 
supplied  to  the  associated  swirlers  and  splash  plate,  since  it  was  thought 
that  these  flows  might  influence  atomization. 

Air  Supply. 

Air  was  supplied  using  a  Metrico  centrifugal  fan  (capable  of  providing 
up  to  1299  CFM  at  a  pressure  of  2  psi) .  The  flow  rate  was  regulated  using  a 
sliding  gate  valve.  Since  air  was  fed  to  the  nozzle  assembly  through  a 
curved  plastic  hose,  honeycomb  metal  flow  straighteners  were  employed 
upstream  of  the  spectacle  plate  inlet. 

The  flow  rate  was  measured  by  means  of  an  orifice  plate,  with  corner 
tappings.  Static  pressure  was  measured  using  a  mercury  manometer  and 
differential  pressure  using  an  inclined  spirit  manometer.  Air  temperature 
was  recorded  using  a  sheathed  chrome 1-alumel  thermocouple  and  an  RS  AD595AD 
thermocouple  amplifier  and  cold  junction  compensator.  The  accuracy  of  flow 
measurement  was  estimated  to  be  ±2%. 

Fuel  Supply. 

Fuel  was  supplied  from  a  steel  pressure  vessel  pressurised  using  a 
nitrogen  cylinder.  The  flow  rate  was  measured  using  a  rotameter  with  an 
accuracy  estimated  to  be  0.005  L/min. 

Fuel  Extraction. 

Spray  droplets  were  collected  in  a  drum  and  fumes  were  extracted  from 
this  by  means  of  a  Coanda  ejector  supplied  with  compressed  air. 

Sizing  Instrument. 

The  sizing  instrument  transmitter  and  detector  were  mounted  on  optical 
trestles  either  side  of  the  fuel  nozzle  and  were  aligned  as  described  in  the 
Malvern  particle  sizer  manual  (3). 

A  300  mm  collecting  lens  was  employed. 


Results. 

Size  distributions  were  measured  through  the  centre  of  the  spray  at 
two  stations  downstream  of  the  injector. 

The  measured  distributions  of  light  energy  were  fitted  with  Rosin 
Rammler  size  distributions.  Results  for  a  range  of  fuel  and  air  flow  rates 
are  shown  in  figures  3  and  4  and  in  table  1. 


Nomenclature. 
a  Radius  of  particle. 

a^  Mean  particle  radius  in  ith  size  range. 

A^  Cross  sectional  area  of  particles  in  the  size  range. 

E  Light  energy  incident  on  particle. 

F  Focal  length  of  Fourier  transform  lens. 

I  Light  intensity  with  particles  in  beam. 

X  Light  intensity  without  particles  in  beam. 

J  Zero  order  Bessel  Function, 

o 

First  order  Bessel  Function. 
i  Optical  path  length, 

m  Number  of  size  ranges. 

N  Number  of  particles. 

N^  Number  of  particles  in  size  range, 

w  Weight  of  particles  radius  a. 

X  Wave  length  of  light. 

&  Particle  density. 
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FIGURE  2  F101  INLET  CALIBRATION  APPARATUS 
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A  STUDY  OF  F101  COUNTER  SWIRLING  INLET. 


Introduction . 

When  rotating  motion  is  imparted  to  a  fluid  upstream  of  an  orifice, 
then  the  flow  leaves  the  orifice  with  tangential  as  well  as  axial  and  radial 
velocity  components.  The  presence  of  this  swirl  sets  up  an  axial  pressure 
gradient  which  causes  reverse  flow  along  the  axis  and  a  recirculation  zone, 
is  formed.  In  a  combustion  chamber  this  is  used  to  stabilise  combustion  in 
the  primary  zone.  The  ability  to  predict  this  phenomenon  will  be  most 
important  for  accurate  combustor  modelling. 

The  object  of  this  work  was  to  make  LDA  measurements  of  velocity  close 
to  the  exit  of  F101  swirlers  and  to  demonstrate  the  ability  of  the  FLUENT 
code  (1)  to  predict  these  results. 

The  F101  Swirlers . 

The  spectacle  plate  and  two  stage  counter  rotating  swirler  assembly 
are  illustrated  in  figure  1.  A  portion  of  the  primary  air  flows  through 
the  swirlers.  These  comprise  a  primary  swirler  of  eight  tangential  slots 
and  a  secondary  guide  van  cascade.  The  remainder  of  the  primary  air 
enters  through  perforations  in  the  structure  to  provide  impingement  cooling 
of  the  splash  plates  and  hot  side  film  cooling  for  the  dome  structure. 

Experimental  Measurements . 

The  calibration  rig  described  in  the  previous  report  (2)  was  again 
used  for  this  work  and  the  LDA  system  is  illustrated  in  figure  2.  The 
latter  employed  an  800  mW  argon  laser  operated  in  the  fundamental  tern 
mode.  Light  is  focused  onto  a  rotating  diffraction  grating  which  serves 
as  a  beam  splitter  and  provides  the  frequency  shifting  required  for  making 
measurements  in  the  highly  turbulent  swirling  flow.  A  mask  blocks  off  all 
but  the  +1  and  -1  order  diffracted  beams  which  are  brought  to  focus  at  the 
measurement  point  in  front  of  the  F101  spectacle  plate.  The  rotating 
grating  is  driven  by  a  servo  motor  and  its  speed  is  monitored  by  means  of  an 
infra  red  transmitter  and  detector. 

Detector  System. 

An  85-200  mm  zoom  lens  was  employed  to  focus  scattered  light  from  the 
measurement  volume  into  a  Malvern  Instruments  photomultiplier  assembly  which 
incorporates  a  pinhole,  interference  filter,  EMI  9502B  photomultiplier, 
amplifier  and  discriminator. 

Clipped  autocorrelation  of  the  signal  was  carried  out  using  a  Malvern 
K7023  digital  correlator  and  correlator  store.  This  is  a  1  bit  clipping 
correlator  and  is  described  fully  in  reference  (3) .  The  device  has  48 
channels  for  storing  the  correlogram.  Sample  times  down  to  0.05  micro 
seconds  can  be  used  and  these  run  consecutively  with  no  dead  time  between. 

The  correlator  was  operated  under  control  from  a  Malvern  MDP7023V 
computing  data  processor.  This  is  essentially  a  microcomputer  running  the 
Malvern  fringes  mode  velocity  and  turbulence  algorithm  (3) .  Data  stored 
in  the  correlator  can  be  down  loaded  to  the  computer  where  it  is  processed 
to  give  values  of  velocity  and  turbulence. 


Velocity  Calculation. 


The  probability  distribution  of  velocity  in  the  turbulent  flow  is 
assumed  Gaussian  with  variance  o.  The  method  employed  to  calculate  the 
mean  velocity,  u,  and  percentage  turbulence 

n  =  (a/u)  x  100%  (1) 

is  that  of  differentiating  the  analytical  expression  for  the  photon 
correlation  function  and  solving  for  the  first  three  turning  points  in  the 
correlogram.  This  yields: 


n  =  V_  i  (R-l)+  — — 


(g2"9l)/(92"g3) 


s/[?(t1+t2+t3)^] 


,  2  1  ,  .  m_ 

1  ~  n  _  222  1  +  2 
m  n  7T 


V  x  V 
o  c 


m  =  (2c)/(a+b) 


a  =  gx  -  FP 


b  =  g2  -  FP 


c  =  g2  "  gl 


where 

m  =  effective  fringe  visibility 

x  =  sample  time 

s  =  fringe  spacing 

n  =  number  of  fringes  in  beam  radius. 

The  meaning  of  t^  and  gi  are  apparent  from  figure  8. 
Gas  velocity  is  then  given  by 


(gas  velocity)  =  Vt  i  C  x  (grating  speed) 
where  C  is  a  constant  dependent  on  the  diffraction  grating  line  spacing. 


Computer  modelling. 

The  FLUENT  code  (1)  was  used  to  model  the  counter  swirling  flow.  The 
boundary  conditions  employed  are  shown  in  tables  1,  2  and  3.  The  finite 
difference  grid  is  shown  in  figure  7(a).  The  inlet  velocity  components  were 
calculated  from  the  combustor  geometry  and  the  known  air  flow  rate  using  the 
swirling  flow  theory  of  Beer  and  Chigier  (4) . 

Modelling  was  carried  out  using  both  the  K  e  model  for  turbulence  and 
the  algebraic  stress  model. 

The  predicted  axial  and  tangential  velocity  components,  at  a  station 
25  mm  downstream  of  the  primary  inlet,  were  compared  with  the  LDA  measurements 

Results. 

Figures  3,  4,  5  and  6  compare  the  predicted  axial  and  tangential 
velocity  components  with  measured  values.  It  can  be  seen  that  using  the 
algebraic  stress  model  for  turbulence  gives  a  considerable  improvement  in  the 
predicted  flow  field. 

Conclusion. 

LDA  measurements  have  been  made  at  the  exit  of  an  F101  counter  swirling 
air  inlet.  The  swirlers  have  also  been  modelled,  using  the  FLUENT  finite 
difference  code,  with  Keand  algebraic  stress  models  for  turbulence.  From 
comparison  of  predicted  and  measured  velocities  it  is  concluded  that  an 
algebraic  stress  model  for  turbulence  gives  a  considerable  improvement  over 
the  Ke  model. 

Nomenclature. 

u  axial  velocity  (m/s) . 

v  radial  velocity  (m/s) . 

w  tangential  velocity  (m/s) . 

WO  denotes  wall  cell. 

I  denotes  inlet. 

O  denotes  outlet. 

S  denotes  symmetry  cell. 

•  denotes  live  cell. 
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THE  F101  COMBUSTOR  TEST  FACILITY  AND  MERCURY  PULSE 
TRACER  MEASUREMENTS. 


Introduction. 

The  finite  difference  programmes  developed  at  Sheffield  for  the 
modelling  of  gas  turbine  combustors,  incorporates  a  stochastic  particle 
tracking  algorithm  to  predict  droplet  trajectories.  Here  individual 
droplets  are  tracked  through  the  combustion  chamber  as  they  experience 
the  predicted  flow  velocities.  At  each  location  a  turbulent  contribution 
is  also  included  whose  magnitude  is  based  on  the  random  choice  of  a  normally 
distributed  variable  simulating  Gaussian  turbulence  and  whose  duration  is 
based  on  an  eddy  decay  model. 

Such  a  procedure  can  also  be  used  for  the  converged  flowfields  with  a 
neutrally  buoyant  fluid  tracer  element.  If  sufficient  elements  are 
tracked  from  a  fixed  input  position  to  the  exit  plane  then  the  resulting 
number  versus  time  is  the  RTD  function  at  the  exit  embodying  the  effects  of 
broadening  due  to  turbulence  and  recirculation. 

By  comparing  these  results  with  measured  RTD  functions  at  the  exit  of 
the  combustor  it  is  possible  to  check  the  accuracy  of  the  predicted  3 
dimensional  flow  field  and  the  particle  tracking  algorithm.  The  results  can 
also  be  used  to  predict  the  blow  off  limits  for  the  combustor,  Ewan  et  al.  (2) . 

This  report  describes  measurement  of  RTD  functions  which  will  be  used 
to  check  the  accuracy  of  F101  combustor  modelling. 

The  F101  combustor  test  facility. 

The  test  combustor  is  a  54°  sector  of  a  full  annular  F101  engine  combustor 
This  comprises  a  combustor  dome,  containing  3  air  blast  fuel  atomizers  and 
associated  inlet  swirlers,  together  with  inner  and  outer  liners  all  of  which 
are  fabricated  from  actual  engine  parts  (figure  2). 

The  chamber  is  mounted  in  a  stainless  steel  vessel.  This  allows 
operation  under  pressure  and  can  be  opened  to  allow  easy  access  for  measurement 
when  operating  at  atmospheric  pressure.  .Air  is  supplied  from  a  four  cylinder 
reciprocating  compressor  at  pressures  of  up  to  100  psi.  High  pressure  and 
realistic  operating  conditions  can  be  achieved  by  discharging  air  from  a 
storage  tank.  Figure  1  shows  the  test  facility. 

Combustor  Instrumentation. 

Total  air  flow  is  measured  by  means  of  a  square  edged  orifice  plate  with 
corner  taps  and  primary  air  flow  by  means  of  a  venturi  meter.  Differential 
and  upstream  pressures  are  measured  using  mercury  manometers.  Air  temperature 
is  monitored  using  chrome 1- a lumel  thermocouples  with  signal  conditioning  by 
means  of  RSAD595AD  amplifier /cold  junction  compensators. 

Mercury  pulse  tracer  method. 

Theory. 

The  residence  time  distribution  (RTD),  denoted  E(t),  is  defined  as  the 
fraction  of  material  at  the  exit  with  a  residence  time  between  t  and  t4dt. 

It  can  be  determined  by  introducing  a  pulse  of  tracer  material  into  the 
chamber  over  a  period  of  time  short  compared  to  the  mean  residence  time  and  by 
monitoring  the  change  of  tracer  concentration  with  time  at  the  exit. 


The  concentration  versus  time  information  can  be  converted  into  the 
form  of  an  RTD  by  normalising.  This  makes  the  result  independent  of  the 
actual  size  of  tracer  pulse.  It  is  usual  to  take  the  average  of  many 
RTD  measurements  to  eliminate  the  effect  of  random  fluctuations  observed 
from  one  experiment  to  the  next.  The  data  in  this  form  is  suitable  for 
comparison  with  computer  model  results. 

Most  tracer  methods  employ  solenoid  operated  valves  to  introduce  steps 
or  pulses  in  a  gas  or  dye  concentration  entering  the  chamber.  However, 
these  methods  often  generate  pulses  with  long  rise  times  compared  to  the  mean 
residence  time.  This  is  a  serious  problem  with  high  intensity  combustors 
where  residual  lifetimes  are  low.  Topps  (1)  has  designed  a  method  in  which 
a  pulse  of  mercury  is  produced  by  passing  an  electric  spark  across  the  ends  of 
two  mercury  amalgam  electrodes.  As  the  mercury  vapour  leaves  the  combustor 
its  concentration  is  monitored.  This  is  done  by  making  use  of  its  ability 
to  absorb  U.V.  light  produced  by  a  low  pressure  mercury  lamp  and  detected  by 
a  photomultiplier.  The  method  has  been  used  successfully  at  Sheffield  on  a 
Lycoming  gas  turbine  combustor  can.  B.C.R.  Ewan  (2).  Tests  have  also 
been  carried  out  in  laminar  pipe  flow  by  means  of  which  it  was  shown  that 
pulses  of  around  1.2  micro  second  duration  could  be  produced. 

Experimental  apparatus. 

The  mercury  spark  electrode. 

This  comprises  a  length  of  twin  bore  refractory  tubing  into  which 
dental  amalgam  has  been  packed.  Details  of  construction  can  be  found  in 
reference  (2) .  The  electrode  was  inserted  into  the  combustion  chamber 
through  a  suitable  hole  drilled  in  the  outer  casing.  The  mercury  vapour 
pulse  was  generated  by  discharging  a  0.1  micro  farad  capacitor  across 
the  tip  after  charging  to  800  V  from  an  EHT  supply. 

Detection  apparatus. 

Tracer  concentration  at  the  exit  was  measured  using  a  low  pressure 
mercury  vapour  lamp  emitting  radiation  at  253.7  nanometres  wave  length  and  an 
EMI  type  9789QB  photomultiplier  tube  with  an  entrance  aperture  of  2  mm 
diameter.  (Hence  a  2  mm  diameter  cylinder  at  the  exit  was  monitored.) 

A  CMB  computer  with  machine  code  program  and  analogue  to  digital 
converter,  were  employed  to  record  the  PM  output  as  mercury  excited  the 
combustor.  Data  acquisition  was  triggered  at  the  instant  of  firing  the 
electrode  using  a  5V  pulse  from  the  capacitor  discharge  circuit.  1024  12  bit 
values  were  recorded  at  0.7  milli-second  intervals.  These  were  then  stored 
on  disc  for  processing  into  RTD  functions  at  a  later  stage. 

Results. 

The  RTD  obtained  from  the  average  of  50  separate  pulse  injections  is 
shwon  in  figure  2.  These  were  obtained  when  the  pulse  was  injected  at  the 
location  of  the  F101  surface  discharge  igniter  and  with  the  combustor 
operating  at  atmospheric  pressure. 

Conclusion. 


Residence  time  distribution  measurements  have  been  made  in  the  FiOl 
gas  turbine  combustion  chamber.  Work  to  obtain  more  data  and  to  simulate 
the  results  numerically  is  continuing. 
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Notation 


*  . 


v  kinematic  viscosity 

u“uj  Reynolds  stress  tensor 

p  pressure 

e  viscous  energy  dissipation  rate 

p  density 

4>..  pressure-strain  tensor 

*  J 

b.  .  anisotropy  tensor 

A  J 

k  turbulent  kinetic  energy 

ci  constant  of  the  linear  part  of  the 

pressure-strain  tensor 

X  constant  of  the  non-linear  part  of  the 

pressure  strain  tensor 

II  second  invariant  of  the  anisotropy  tensor 

III  third  invariant  of  the  anisotropy  tensor 

c£2  constant  of  the  model  for  energy  dissipation 
TI  mean  velocity  in  the  xi  direction 

•t  normalised  time 

t  elapsed  time 

ai  constant  equal  to  ci~1 

u2  streamwise  component  of  the  Reynolds 

stress  tensor 

v2  lateral  component  of  the  Reynolds 


stress  tensor 


'component  of  the  Reynolds  stress  tensor 
% 

ifc  the  x3  direction 
perturbation  parameter 

♦ 

parameter  that  measures  the  departure 
from  axi symmetry 
production  tensor 

half  the  trace  of  the  production  tensor 

principal  components  of  the  Reynolds 

stress  tensor 

space  coordinates 

grid  Reynolds  number  ^ 

principal  component  of  the  Reynolds 

stress  tensor 

principal  component  of  the  Reynolds 
stress  tensor 

principal  component  of  the  Reynolds 
stress  tensor 


superscript 
average  value 


right  hand  side  are  the  pressure-strain  and  the  dissipation  terms 
respectively.  The  expression  for  the  pressure  strain  term  follows 
from  the  fact  that  in  homogeneous  turbulence 


Although,  without  analytical  proof,  it  can  be  said  on  physical 
grounds,  that  the  role  of  the  pressure-strain  term  is  to  force  an 
equilibration  of  the  diagonal  components  of  the  Reynolds  stress 
tensor,  nothing  can  be  concluded  from  equation  (l)  about  the 
direction  of  the  energy  transfer  and  the  speed  at  which  the 
tendency  towards  equipartition  is  evolving.  The  dissipation  term 
in  the  above  equation  is  usually  replaced  by  its  isotropic  part 
with 

e  =  v 

L_  _j 

The  transport  equation  for  the  viscous  dissipation  rate  in 
unstrained,  high  Reynolds  number,  homogeneous  turbulence, 
neglecting  diffusional  effects  of  pressure  fluctuations  can  be 
expressed  as: 


au^^  j 


where,  following  Hanjalic  and  Launder  (1972),  the  first  term  on 
the  right  hand  side  represent  the  generation  rate  of  vcrticity 


fluctuations  through  the  self -stretching  action  of  turbulence, 
and  the  second  term  is  the  decay  of  the  dissipation  rate  through 
the  action  of  viscosity. 

Unfortunately  there  is  no  a  direct  solution  for  these 
complex  equations  even  in  the  case  of  homogeneous  turbulence  and 
therefore  only  qualitative  information  can  be  inferred  from 
these.  To  advance  in  the  solution  of  the  problem  it  is  necessary 
to  introduce  modelling  assumptions. 

For  the  transport  equation  of  the  Reynolds  stress 
components  the  explicit  appearance  of  the  pressure-strain  term 
satisfies  a  quasi-Poisson  equation  as  shown  by  Chou  (1945).  The 
solution  for  this  Poisson  equation,  for  the  case  of  unstrained 
flows  shows  that  the  pressure-strain  originates  from  the 
interaction  between  turbulence  components.  On  physical  grounds 
Rotta  (1951)  suggested  that  this  term  may  be,  as  a  first 
approximation,  proportional  to  the  anisotropy  tensor  b. ^ 
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where  cj  is  a  constant  of  proportionality,  and  k  =  ^u.u.  the 


turbulent  kinetic  energy. 


Using  different  approaches  to  the  modelling  of  this 
term,  lumley  (1977)  Speziale  (1980),  Launder  (private 
communication)  have  expanded  this  expression  incorporating  non¬ 
linear  terms.  On  the  whole,  based  on  the  current  pressure-strain 
models  proposed,  equation  (1)  can  be  expressed  in  the  following 
general  form: 

DuTuT 

Dt  °  =  "°iebij  +  Xe^bikbkj  “  3II6ij)  ”  3E6ij  ^ 

where  II  =  b^b^  is  the  second  invariant  of  the  anisotropy 

tensor,  cj  and  X  can  be  either  constants  or  functions  of  some 
relevant  parameters . 

For  the  equation  of  viscous  dissipation  rate  Rodi  • 

(1971)  has  argued  that  both  terms  in  the  right  hand  side  of 

equation  .,(2)  should  be  modelled  together.  These  terms  are 
governed  by  the  dynamics  of  the  transfer  of  energy  from  low  to 
high  wave  numbers  and  so  they  are  independent  of  viscosity;  on 
dimensional  considerations  it  can  be  suggested  that: 


du.au. bu,  |  a2u. 

> — = - = - +  v - — 

dXjdxk5xk  j  dx^dXj 
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t2k 


where  c^  is  a  constant.  Lumley  and  Khajeh-Nouri  (1972)  have 
suggested  independently  that  the  factor  c g  should  be  a 
dimensionless  invariant  function,  and  that  this  invariant 


function,  in  the  general  case ,• contains  the  entire  mechanism  of 
production  and  destruction  of  dissipation. 


Thus,  based  on  the  currents  models,  the  equation  for 
the  viscous  dissipation  rate  is 


-  _c  £_ 

Dt  “  e2k 


(4) 


Equations  (3)  and  (4)  constitute  a  closed  system  and 
can  he  solved  using  numerical  procedures.  The  constants  involved 
in  the  modelling  are  tuned  in  the  light  of  experimental  data.  But 
since  the  modelling  of  the  rate  of  energy  dissipation  equation  is 
usually  based  on  simplifying  assumptions,  the  conclusions  can  be 
sometimes  misleading.  A  rather  different  approach  is  adopted  here 
to  obtain .the  variation  of  the  components  of  the  Reynolds  stress 
tensor,  along  the  flow. 


2 .-  Equations  for  the  Invariants . 

It  may  be  shown  empirically  that  the  decay  of  turbulent 
kinetic,  energy  implies,  for  an  initially  anisotropic  turbulent 
flow  without  mean  velocity  gradients,  a  decay  in  the  state  -of 
anisotropy  and  consequently  a  decay  oi  the  second  invariant  of 
the  anisotropy  tensor  (this  second  invariant  can  be  regarded  as  a 
measure  of  the  anisotropy  of  the  flow),  and  also  a  decay  of  the 
absolute  value  of  the  third  invariant.  By  constructing 
differential  transport  equations  for  these  invariants  (see 
Appendix  A)  it  would  be  possible  to  obtain  more  information  of 
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To  solve  these  equations  the  normalised  time  t  is  taken  from  the 
experiments.  In  this  way  the  problem  of  solving  a  tensor  equation 
with  three  components  and  a  scalar  equation  for  the  energy 
dissipation  with  the  elapsed  time  as  the  independent  variable  has 
been  -transformed  into  the  solution  of  three  scalar  equations  with 
the  normalised  time,  which  is  related  to  the  natural  time  of 
decay,  as  the  independent  variable.  Thus  by  solving  the  equations 
for  the  turbulent  kinetic  energy,  the  second  and  third  invariant 
of  the  anisotropy  tensor,  the  implications  of  the  assumptions 
involved  in  the  modelling  of  the  energy  dissipation  has  been 
eliminated  and  the  effects  of  the  pressure-strain  term  isolated. 
The  physical  meaning  of*  the  second  invariant  II  is  that  it  is  a 
measure  of  the  state  of  anisotropy  of  the  flow,  and  the  role  and 
physical  meaning  of  the  third  invariant  in  the  return  to  isotropy 
process  shall  be  dealt  in  the  next  section. 


^5.-  Cubic  equation  for  the  Principal  Reynolds  stress  components . 

For  a  given  variation  of  the  invariants  with  respect  to 
the  normalised  time,  the  principal  components  of  the  Reynolds 
stress  tensor  can  be  recovered  from  the  following  algebraic 


relations: 


\i  - 0 


bikbki  *  11 


bikbkjbji  =  111 


For  the  axisymmetric  case  the  diagonal  components  of  the 


t>22  =  1>3  3  =  lb  ,  with  b>0  ,  giving  II  =  ^b2  and  III  =  ~b?  so 
that  III<0  or  as  bn  =  b  ,  b2  2  =  ^>3  3  =  -*^b  ,  with  b>0,  giving 

II  =  ^b2  and  III  =  ^b3  thus  III>0,  it  can  be  seen  from  here  that 

2  4 

III  is  positive  when  two  components  of  the  Reynolds  stress  tensor 

2 

are  less  than  the  isotropic  level  and  III  is  negative  when  two 

3 

components  are  larger  than  the  isotropic  level  (a  similar 
explanation  was  first  given  by  Gence  and  Mathieu  (1980),  the 
third  invariant  is  positive  when  the  component  of  uhe  Reynolds 
stress  tensor  along  the  symmetry  axis  is  greater  than  the  two 
others,  and  is  negative  in  the  opposite  case).  The  solutions  of 
the  algebraic  relations  for  this  particular  case  are: 
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where  the  sign  to  be  taken  follows  frojn  the  above  discussion, 
depending  on  the  sign  of  the  third  invariant  III. 

Moreover  for  the  general  case  when  three  components  of 
the  Reynolds  stress  tensor  are  distinct  a  cubic  equation  may  be 
obtained: 

x3  -  illx  -  illl  =  0  (15) 

with  x  representing  the  individual  components  of  the  anisotropy 
Tensor.  For  the  problem  to  have  a  physical  meaning  the  three 
roots  of  equation  (15)  must  be  real  since  each  root  corresponds 
to  each  diagonal  component  of  the  anisotropy  tensor.  From  the 
algebraic  solution  of  the  cubic  equation 

F(x)  =  ax3  +  3bx2  +  3cx  +  d  =0 
with  Q  =  ac  -  b2  and  R  =  l(3abc  -  da2)  -  b3  it  follows  that 

(i)  If  -Q3  +-R 2  -  Q  ie.  6III2  -  II3  =  0  ,F(x)  has 
three  real  roots  of  which  at  least  two  are  equal,  this  is  the 
axisymmetric  case  and  a  strong  realizability  condition  for 
axisymmetry  appears  in  terms  of  the  invariants,  ie.  6III2  - 
II3  =  0 


-  9  - 


Us**®* 


(ii)  If  Q3  +  R2  =0  ie.  6III2  -  II3  >  0  F(x)  has  one 
real  root  only  and  the  other  two  are  complex,  this  solution  has 
no  physical  meaning. 

(iii)  If  Q3  +  R2  =0  ie.  6III2  -  II3  <  0  F(x)  has 
tfiree  distincts  real  roots  which  correspond  to  the  three 
principal  velocity  correlations. 

The  solution  for  the  principal  components  of  the 
Reynolds  stress  tensor  is: 


=  kli  + 


_  i  II! 
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(16) 


where  n  =  0,  2,  4  ;  when  III>0  the  upper  sign  is  used  and  when 
III<0  the  lower  sign  is  used  (the  sign  follows  from  the  sign  of  R 
in  the  general  solution  for  the  cubic  equation,  which  in  our  case 
turns  out  to  be  the  sign  of  the  third  invariant) 

Observing  these  simple  analytical  solutions  for  the 
energy  components  a  "physical  picture  of  the  evolution  of  these  in 
return  to  isotropy  can  be  suggested.  It  can  be  noted  that  the 
sign  of  the  third  invariant  remains  constant  during  the  process 
of  return  to  isotropy.  A  change  of  sign  means  that  the  turbulence 
have  passed  through  a  state  where  there  are  only  two  equal  and  of 
opposite  sign  nonzero  components  for  the  anisotropy  tensor  so 
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that  111=0,  (the  diagonal  components  of  the  anisotropy  tensor  in 
this  state  are  of  the  form  tu  =  b  ,  b22  =  “b  »  and  b33  =  0) . 
This  phenomenon  is  not  physically  probable  without  the  external 
production  of  energy  through  the  action  of  velocity  gradients  and 
it  follows  therefore  that  the  variation  of  the  angles 


U_-1  [il3  -  6III 2] 2  . 
3*  !  6m2  ! 


with  n  =  0,  2,  4  cannot  be  equal  or  more  than  30°  ,  therefore  no 

change  of  the  sign  of  the  cosine  of  these  angles  is  expected 

during  the  process.  When  the  third  invariant  is  negative  the 

smallest  component  has  the  largest  absolute  value  of  the 

respective  cosine  term  in  the  solution  of  the  cubic  equation,  and 

2 

it  is  approaching  isotropy  from  below  the  isotropic  level  — k. 
This  component  is  gaining  energy  from  the  other  two  larger 
components,  (both  are  approaching  isotropy  from  above  the 
isotropic  level)  in  a  proportion  which  may  be  related  to  the 
absolute  value  of  their  respective  cosine  terms.  Therefore  it  is 
likely  that  there  will  be  an  increase  of  energy  in  the  smaller 
component  and  the  process  of  return  to  isotropy  will  be  faster. 
When  the  third  invariant  is  positive  the  picture  is  quite  the 
contrary.  The  largest  component  has  now  the  largest  absolute 
value  of  the  cosine  term  and  it  is  giving  energy  to  the  other  two 
smaller  components  (which  are  approaching  isotropy  from  below  the 
isotropic  level).  If  the  relative  amount  of  energy  shared  by  the 


1 1 


smaller  components  is  taken  to  be  proportional  to  the  absolute 
value  of  their  respective  cosine  terms,  the  smaller  the  component 
the  greater  is  the  amount  of  relative  energy  accumulated  from  the 
largest  component.  But,  because  of  the  sharing  it  is  unlikely  for 
the  smallest  component  to  have  a  significant  increase  of  energy 
during  the  process  and  the  return  to  isotropy  is  likely  to  be 
slower  than  when  the  third  invariant  is  negative.  With  this 
mechanism  in  mind  it  is  seen  that  the  direction  of  the  flow  of 
energy  and  possibly  the  relative  amount  of  energy  transfered 
among  the  components  is  governed  by  the  sign  of  the  third 
invariant,  the  sign  of  the  cosine  terms  in  the  solution  of  the 
cubic  equation  and  the  absolute  value  of  their  respective  cosine 
terms.  It  should  be  mentioned  that  the  influence  of  the  third 
invariant  in  the  dynamics  of  return  to  isotropy  was  suggested  by 
the  Genes  and  Mathieu  experiment  (1980). 

Another  observation  from  the  solution  of  (15)  is  that 
there  appears  to  be  a  realizability  condition  for  the  existence 
of  unstrained  homogeneous  turbulence  in  terms  of  the  invariants 
which  can  be  expressed  as  6III2  -  II 0;  but  since  in  theory 
any  turbulent  flow  may  be  characterized  by  its  invariants  this 
general  realizability  condition  should  be  satisfied  by  any 
proposed  model  of  turbulence.  The  simplicity  of  this  expression 
constitutes  a  very  simple  test  for  the  phenomenological  models. 

4.-  Some  solution  for  the  equations  of  the  Invariants. 
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.  The  next  step  ..  in  this  analysis  is  to  solve  the 
equations  for  the  invariants  as  a  function  of  the  normalised  time 


4.1-  Neglecting  non-linear  effects  (  X  =  0  ) 

When  ci  is  a  constant  we  have  Rotta's  proposal  (1$51), 
and  the  analytical  :  ..-tions  for  II  and  III  are 


II  =  line 


-2(  Cl  -  1)' 


(17) 


III  =  III0e 


-3(  cj  -  1)‘ 


(18) 


with  II0  and  III0  being  the  initial  values  for  the  second  and 


II3 


third  invariants.  From  these  solutions  we  obtain  jffT  =  constant. 


which  is  only  true  for  homogeneous  axisymmetric  turbulent  flows 
(constant=6)  or  at  most  when  the  departure  from  axisymmetry  is 
small.  Moreover  ci  can  be  estimated  using  the  experimental  data: 


ln(  S-) 


c  i  =  1  + 


1 


II, 


21"<  ft 


with  k<k0  ,  II<II0  everywhere,  and  since  decays  means  that  the 
values  of  the  second  invariant  and  the  kinetic  energy  in  the 
process  are  less  than  their  respective  initial  values  it  follows 
that  cj>  1  .0.  (this  can  also  be  seen  from  the  solutions  (17)  and 
(18)  where  if  ci<1 .0  there  would  be  an  increase  in  the  anisotropy 
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of  the  flow  without  the  action  of  a  generation  term) 

* 

When  one  uses  the  data  of  Uberoi  (195  )  (anisotropic 
decaying  turbulence  submitted  to  axisymmetric  plane  strain  with  a 
contraction  ratio  of  4:1)  to  estimate  the  constant  cj  ,  it 
becomes  apparent  that  there  is  not  a  unique  value  which  gives  the 
same  agreement  for  the  three  different  contraction  ratios  of  the 
axisymmetric  duct  used  in  the  experiments.  Furthermore,  the  value 
of  Ci  varies  along  the  development  of  the  flow:  1 .2  ^  ci  2.3. 
As  a  first  approximation  ci  can  be  taken  a  constant  with  ci  =2.0 
giving  reasonably  agreement  within  the  uncertainty  of  the 
experimental  data.  Figures  1,  2,  and  3  shows  this  satisfactory 
agreement. 


When  X  =  0  and  ci  is  not  a  constant.  For  this 
situation  equations  (7)  and  (8)  may  be  reduced  to: 

(3III)~i  =  (211)jp  (19) 

and  solving,  one  obtains 
II3 

Yffj  ~  constant 

Again  the  same  condition  as  for  cj  being  a  constant  is  obtained. 
Therefore  any  improvement  in  the  modelling  of  the  pressure-strain 
term,  for  the  return  to  isotropy,  in  expressing  cj  as  a  function 
of  some  relevant  parameters  (Lumley  (1978))  is  still 
unsatisfactory  to  model  non-axisymmetric  flows. 


4.2  Non-linear  term  in  the  pressure-strain  tensor 


Focusing  attention  on  the  complete  expression  for  the 
pressure-strain  tensor  in  the  return  to  isotropy  process  and  for 
simplicity  taking  cj  and  A  to  he  constants,  the  constraints 
imposed  by  the  linear  model  can  be  avoided. 

(i)  Axisymmetric  flows 

A  straightforward  analytical  solution  is  obtained  in 
the  axisymmetric  case: 


when  substituted  into  equation  (8)  one  obtains 
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dx 
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which  can  be  easily  integrated  to  give 


II  = 


II0(  c,-1 )2e 


-2(  c  j—1  )  x 
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and  similarly  for  III 
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vhen  III<0  the  upper  sign  is  used  and  vice  versa 

Again  by  taking  Uheroi  (1956)  experiment,  it  is 

possible  to  estimate  the  value  of  second  constant  A  by  fixing  cj . 
Along  the  flow  the  variation  of  A  is  small  for  a  fixed  ci  thus 

||  1  f 

A  J  II0  is  plotted 

against  (  cj-1 )  a  straight  line  relationship  between  the 

constants  is  obtained  (see  figure  4).  It  is  noted  here  that  the 
sign  of  A  depends  on  the  sign  of  the  third  invariant.  Figures  5, 

6,  and  7  show  the  prediction  for  the  Uberoi  (1956)  experiment 

with  different  values  for  ci  and  A  obtained  from  Figure  4.  An 
improvement  with  respect  to  the  linear  model  can  be  noted. 


( ii )  Non-axisymmetric  flows . 

Vhen  the  flow  is  not  axisymmetric  an  approximate 
analytical  solution  can  be  obtained  using  perturbation  methods. 
It  can  be  assumed  that  A  is  small  parameter  and  that  the  solution 
for  the  invariants  is  also  a  function  of  this  parameter  together 
with  the  constant  cj. 

Writing  equations  (10)  and  (11)  in  the  form: 
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ai  =  (  cj-1),  where  IIq  and  IIIo  are  the  initial  values  for  the 
invariants  at  t  =  0  and  rearranging  one  obtains 


u. 


fl  +  5&1II  +  6af II  -  X2II2  =  0 

Normalising  the  invariants  with  respect  to  their  initial  values 
and  putting  w2  =  X2IIq  .with  0  <  w  <  1  the  perturbation  parameter 


so  that  0  <  X  < 


1 


Jn7 


equation  (22)  becomes: 


II  +  II  +  6a?II  -  u>2II  2  =0 
n  x  n  1  n  n 


t^is  equation  the  following  transformation  is  applied: 
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Equation  (23)  now  reduces  to 


1  “  ialT 

X  -  lai2X  -  u)2e  2  X2  =  0 


which  can  be  easily  solved  using  perturbation  methods,  (see 
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Finally  when  the  initial  conditions  are  introduced  the 
solutions  for  the  invariants  are: 


x2Ilg“]  -4a!  t 
-2ZT\e 


(24) 


r  2xllg  -3a! x 
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Expressions  (24)  and  (25)  indicate  that  X  and  ai  are 
important  parameters  which  influence  the  speed  at  which  the 
return  to  isotropy  is  taking  place.  When  the  invariants  were 
normalised  using  the  initial  conditions,  X  is  related  to  the 


perturbation  parameter  w  as  X  =  w/  J  II  o.  for  0  .<  u  <  1  .  Then  Once 
w  is  fixed  the  value  of  X  is  obtained  from  the  initial  condition 
of  the  second  invariant.  In  section  4  it  has  been --remarked,  that 
the  third  invariant  has  an  important  influence  on  the  transfer  of 
energy  among  the  components  and  also  on  the  speed  at  which  the 
return  to  isotropy  takes  place;  it  is  then  possible  to  reflect 
how  fast  or  slow  the  process  of  equilibration  of  energy  is  taking 
place  through  the  constant  aj.  For  axisymmetry,  in  the  linear 
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model,  ax  -1.0  gives  satisfactory  agreement  with  the  data  and  the 

perturbation  solution' is  around  this  state. 

A  new  parameter  which  measures  the  initial  degree  of 

departure  from  axisymmetry  can  be  defined  as: 

6IIlg 

--uT 

so  that  c  =  0  for  axisymmetry  and  Q— >  1  as  the  departure  from 
axisymmetry  increases  On  an  empirical  basis  the  following 
expression  for  ai  can  be  put  forward: 

—1 OqIII o  (26) 

ai  =  e 

-fom0  ^  . 

(or  more  generally  =  e  w'th  f  being  a  factor  of 

proportionality  ),  thus 

j~  <  1  for  III0  >  0 

ai  =  |  1  for  axisymmetry 
j  >  1  for  Illo  >  0 
i_ 

5.-  Results  for  non-axisymmetri c  flows 

The  above  analysis  is  now  applied  to  three  experiments 
on  the  return  to  isotropy:  Tucker  and  Reynolds  (1968),  Gence  and 
Mathieu  (1980),  Le  Penven,  Gence  and  Comte-Bellot  (1984).  Table  1 
presents  the  conditions  used  in  the  calculation  of  these 


experiments . 


5 « 1  Tucker  and  Reynolds  Experiment  ( 1 968) 

In  this  experiment'  homogeneous  grid-generated 
turbulence  was  subjected  to  a  uniform  plane  strain  through  the 
action  of  a  distorting  duct  with  a  rectangular  constant  cross- 
sectional  area,  measurements  which  reveal  the  tendency  of  the 
flow  to  become  less  anisotropic  were  taken  in  the  parallel  duct. 
The  sign  of  the  third  invariant  is  negative  and  the  component 
with  less  energy  is  increasing  along  the  process  of  return  to 
isotropy.  Le  Penven  et  al  (1984)  pointed  out  that  this  is  a 
visible  effect  of  the  influence  of  the  third  invariant  in  the 
dynamics  of  equilibration  of  the  energy  components.  Figures  8(a) 
and  8(b)  show  the  comparison  between  experimental  data  and  the 
analytical  solutions  (21)  and  (22)  together  with  that  of  the 
linear  model  (equation  18).  Good  agreement  is  obtained  with  the 
perturbation  solution  in  particular  for  the  smallest  component 
which  increases  its  energy  along  the  parallel  duct,  demonstrating 
the  superiority  of  the  non-linear  model.  It  is  necessary  to 
mention  that  the  total  decay  time  was  taken  from  the  experiment 
but  the  decay  time  at  each  measurement  in  the  flow  was  estimated 
from  the  least  square  fit  to  the  experimental  data,  since  the 
experimental  turbulent  kinetic  energy  does  not  give  a  smooth 
decay  due  to  scatter. 

5 ‘2  Gence  and  Mathieu  Experiment  ( 1 980) 

These  authors  have  considered  the  return  to  isotropy 


> T 


phenomenon  after  the  homogeneous  turbulent  flow  has  been 

submitted  to  two  successive  plane  strains  whose  principal  axed 

are  shifted  by  an  angle  a  .  Two  important  conclusions  are 

expressed  in  their  paper,  (a)  The  orientation  of  the  principal 

axes  of  the  Reynolds  stress  tensor  does  not  change  during  the 

return  to  isotropy  process  and  is  independent  of  the  history  of 

the  turbulence  previous  to  the  starting  point  of  the  return  to 

isotropy,  and  furthermore,  the  axes  of  the  pressure  strain  term 

and  the  axes  of  the  Reynolds  stress  tensor  are  aligned  during  the 

whole  process,  and  (b)  the  third  invariant  is  an  important 

parameter  which  influences  the  mechanics  of  the  return  to 

isotropy,  the  process  being  slower  or  faster  according  to  the 

sign  of  the  third  invariant.  Pig  9(a)  to  11(b)  present  the 

comparison  between  the  linear  model  (equations  (17)  and  (18)), the 

perturbation  solution  and  the  experimental  data.  The  fact  that 

there  is  no  change  of  orientation  of  the  principal  axes  of  the 

Reynolds  stress  tensor,  when  the  angle  between  both  strains  were 

shifted,  was  utilized  in  the  calculations,  and  the  principal 

components  were  calculated  from  the  solution  of  the  equations  for 

the  invariants.  For  the  angles  a  =  0  and  a  =  ~  excellent 

4 

agreement  is  obtained  with  the  perturbation  solution  (  being 

superior  to  the  linear  model).  When  the  angle  o  =  3^  if  seems 

4 

that  the  nonlinear  solution  does  not  reproduce  adequately  the 
above  agreement,  but  this  can  be  attributed  to  the  fact  that  the 
initial  conditions  which  determine  the  values  of  the  constants  aj 
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and  X  and  these  do  not  reflect  the  overall  state  of  the  flow  (see 
Table  1).  The  data  of  Genq^  and  Mathieu  (1980)  clearly  shows  that 
the  second  invariant  for  o  =  3j  does  not  decay  smoothly  in  the 
first  three  measurement  locations  of  the  non-distorting  duct. 
Provided  that  the  experimental  data  does  not  reflect  major  errors 
a  possible  explanation  for  this  situation  is  seen  when  the 
balance  of  energy  of  the  smaller  component  can  be  visualised;  the 
smaller  component  is  gaining  energy  from  the  largest  and  losing 
some  energy  by  the  dissipative  action  of  viscosity,  but  just 
after  the  strain  is  suppressed  the  gain  of  the  smallest  component 
is  compensating  the  decay  and  therefore  the  anisotropy  of  the 
flow  is  not  decaying.  The  initial  time  of  decay  should  then  be 
displaced  to  the  point  when  there  is  a  clear  decay  of  the 
anisotropy.  Figure  11.1(a)  and  11.1(b)  shows  this  correction  and 
the  agreement  iB  excellent. 

5. 3  Le  Penven ,  Gence  and  Comte-Bellot  Experiment  (1984) 

This  experiment  is  interesting  in  the  confirmation  of 
the  ideas  put  by  Gence  and  Mathieu  (1980).  In  this  experiment  the 
authors  have  observed  the  evolution  of  the  kinetic  energy* .second 
invariant  and  the  Reynolds  stress  components,  under  the  influence 
of  the  third  invariant  with  different  sign.  Conclusions  obtained 
are  similar  to  the  Gence  and  Mathieu  (1980)  experiment  but  they 
go  further  in  giving  a  physical  interpretation  of  the  role  of  the 
third  invariant  III,  and  emphasize  its  importance  in  the 


modelling  of  turbulent  flows. 

Figures  12(a)  to  13(b)  present  the  predictions  for  this 
experiment,  when  the  third  invariant  is  negative  good  agreement 
is  obtained  with  the  perturbation  solution.  The  linear  solution 
appears  to  be  inferior  in  the  prediction  of  the  smallest 
component.  When  the  third  invariant  is  negative  both  solutions 
seems  to  be  unsatisfactory,  but  again  using  the  same  argument  as 
for  the  third  case  in  the  experiment  of  Gence  and  Mathieu  (1980) 
and  shifting  the  initial  time  of  decay,  the  predictions  agree 
with  the  experiment  with  the  non-linear  solution  being  better 
than  the  linear  one. 

6.-  Conclusions 

The  present  analysis  shows  that  the  discrepancies 
between  experiment  and  prediction  with  the  current  models  of 
turbulence,  even  for  the  simpler  case  of  return  to  isotropy,  can 
be  largely  due  to  the  modelling  of  the  rate  of  dissipation  of 
energy  equation.  In  the  return  to  isotropy  process  the  effects 
from  the  modelling  of  the  pressure-strain  tensor  and  the  e- 
dissipation  may  be  isolated  by  using  a  normalised  time  related  to 
the  natural  time  of  decay. 

Rotta’s  proposal  (1951)  and  those  proposals  where  the 
factor  of  proportionality  in  the  modelling  of  the  pressure-strain 
term  is  not  a  constant,  are  still  unsatisfactory  to  model  non- 
axisymmetry  flows.  Nevertheless,  the  linear  model  with  cj  being  a 


constant  gives  reasonable  agreement  and  the  predictions  can  be 
improved  by  introducing  a  non-linear  term  with  a  second  constant 
X.  In  this  way  it  is  expected  to  avoid  the  constraint  of 
axisymmetry  imposed  by  linear  modelling 

A  physical  picture  of  the  transfer  of  energy  among  the 
components  of  the  Reynolds  stress  tensor  is  suggested  on  the 
basis  of  the  solution  of  the  cubic  equation  for  the  principal 
components  of  the  Reynolds  stress  tensor;  when  the  third 
invariant  is  negative  the  two  larger  components  are  giving  energy 
to  the  smaller  consequently  the  smallest  component  tend  to 
increase  its  energy;  while  when  the  third  invariant  is  positive 
thf'  two  smaller  component  are  gaining  energy  from  the  largest  and 
the  return  to  isotropy  is  slower  than  the  former  case. 
Furthermore  Le  Penven  et  al.  (1984)  point  out  that  when  the  third 
invariant  is  positive  there  may  be  a  return  to  axisymmetry,  but 
with  the  suggestions  of  section  4»  it  is  possible  to  explain  that 
this  phenomenon  is  not  likely  to  happen  in  all  cases  and  depend 
on  the  particular  initial  configuration  rather  than  the  sign  of 
the  third  invariant;  if  the  two  smaller  (greater)  initial  energy, 
component  are  nearly  equal,  there  may  be  a  return  to  axisymmetry 
from  "below"  ("above")  the  isotropic  level  when  the  third 
invariant  is  positive  (negative),  due  to  the  direction  and  the 
relative  amount  of  the  flow  of  energy  transfered. 

Finally  it  should  be  emphasized  that  the  six  cases  of 
these  three  experiments  were  well  reproduced  by  the  perturbation 


References 


1 Chou,  P.  Y.  (1945)  "On  velocity  correlations  and  the  solution 
of  the  equations  of  turbulence  fluctuation".  Quart.  Appl.  Math.  , 
3,  38-54. 


2.-  Gence,  J.  N.,  and  Mathieu  J.»  (1980)  "The  Return  to  isotropy 
of  an  homogeneous  turbulence  having  been  submitted  to  two 
successive  plane  strains",  J.  Fluid  Mech.,  Vol  101,  555-566 


3*-  Le  Penven,  L.,  Gence,  J.N.  and  Comte-Bellot ,  "On  the  Approach 
to  isotropy  of  Homogeneous  Turbulence:  Effect  of  the  Partition  of 
the  Kinetic  Energy  among  the  Velocity  Components",  (1984) 
Fundamentals  of  Fluid  Mechanics  Conference. 


4.-  Hanjalic,  K.  &  Launder,  B.  E.  "A  Reynolds  stress  model  of 
turbulence  and  its  application  to  thin  shear  layer  flows".  J. 
Fluids  Mech.  52,  609-638. 


5. -Lumley,  J.  1. ,  "Computational  Modelling  of  Turbulent  Flows", 
Advances  in  Applied  Mechanics,  Vo.  18.  (1978) 


6.-  Lumley,  J.  L. ,  and  Newman,  G.  R.,  "The  return  to  isotropy  of 
homogeneous  turbulence",  J.  Fluid  Mech.  (1977),  Vol  82,  pp.  161  — 
178 


7.-  Lumley,  J.  L. ,  &  Khajeh-Nouri ,  B.  (197?)  "Computational 
Modeling  of  Turbulent  Transport",  Adv.  in  Geophysics.  A-18,  169- 
192. 


8.-  Rodi,  W.  (1971)  "On  the  equation  governing  the  rate  of 
turbulent  energy  dissipation".  Mech.  Engng.  Dept.  Imperial 
College,  Rep.  TM/TN/A/14 


9.-  Rotta.  J.  C.  (1951)  "Statistische  Theorie  nichthomogener 
Turbulenz.  Z.  Fhys.,  129,  547-572. 


10.-  Tucker,  H.  J.  ,  &  A.  J.  Reynolds  (1968).  "  The  distortion  of 
turbulence  by  irrotational  plane  strain",  J.  Fluid  Mech.,  3?,  657 


-  26  - 


APPENDIX  A 


Equation  for  the  invariants 


The  equations  for  the  invariant  of  the  anisotropy 
tensor  are  obtained  from  the  transport  equation  of  the  Reynolds 
stress  tensor  and  from  the  transport  equation  of  the  turbulent 
kinetic  energy 


Du  u  2 

Dt  *ij  Yij  3  ij 


h  =  p-e  (2) 

where  P^  ,  Y^  ,  and  e  are  the  production  tensor,  the  pressure- 
strain  tensor  and  the  viscous  dissipation  of  energy  respectively. 
The  production  P  and  the  turbulent  kinetic  energy  k  are  half  of 
the  trace  of  the  production  and  the  Reynolds  stress  tensor 
respectively 


The  anisotropy  tensor  and  its  invariant  are  defined  as 


follows 


UiUj 

k 


2, 

36i0 


11  =  bik\i  and  111 


bikbkjb;ji 


then 


+  lb,  ,)) 


=  P.  .  +  Y.  .  -  •=-e6_.  , 


,  r'.:.v.wv.v  ,\ 


and  using  equation  (2) 


Db, 


(  bii  +  f6lj)(p  -  «>  +  kHU  “  p"  'i 


ij  V1!  '  ’3e6D 


Db . 

2  bijbij(  p  "  E)  +  2kbijDt  =  2bijPij  +  2bi3fij 


it  thus  follows 


klr  ■  2(  bijPij  +  t13Ii3  )  -  2II<  P  -  «  >  (5) 

In  a  similar  way  the  equation  for  the  third  invariant  can  be 
obtained 

Hsr1  •  5(  bikbk/u  +  bikbk3Ti3  >  -  3III<P  - « >  - 2IIP  <4> 


In  the  return  to  isotropy  process  P^  =  0,  and  once  inserting  the 
expression  for  the  pressure-strain  tensor. 


Yij  =  "e^  Clbij  “  ^  bikbkj  “  5Il6ij^» 
the  equations  (3)  and  (4)  are  reduced  to 

=  -2|  (  cj  -1)11  -  XIII) 

§F  =  (  Cl  -Din  -  in"  ) 


(5) 

(6) 
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APPENDIX  B 


Solution  for  the  Invariants  using  Perturbation  Methods 


Equations  (10)  and  (11)  of  section  3  can  be  expressed 
in  a  simpler  form  as: 


x  =  -2aix  +  2\y 


y  *  -3aiy  +  i\x2  (2) 

where  x=III ,  y.=III  aj  =  Cj  -1  ,  A  and  ci  being  the  constants 
connected  to  the  non-linear  and  linear  part  of  the  pressure- 
strain  term  respectively.  The  differentiation  is  with  respect  to 
the  defined  normalised  time. 

Taking  derivative  with  respect  to  the  noemalised  time 
again  and  rearranging  the  equations  (1)  and  (2)  one  can  obtain 

x  +  5ajx  +  6ajx  -  X2x2  =  0  (3) 

Normalising  the  second  invariant  with  respect  to  their  initial 
values,  X  =  x/x  equation  (3)  becomes: 

°  9 

X  +  5aiX  +  6a?X  -  X2x0X2  =0 

Now  defining  u  as  the  perturbation  parameter  so  that  w2  =  A2xQ 
for  0<  u  <1  then 


X  +  5ajX  +  6a?X  -  w2X2  =  0  (4) 

To  solve  this  equation  the  following  transformation  is 


applied: 


X  =U  e  z 


equation  (4)  reduces  to 

_5aix 

U  -  iafU  -  w2e  2  U2  =  0. 

Assuming  that  the  solution  may  be  expressed  in  the  form 
TJ  =  TJq  +  ojOj  +  u2U2  +«.. 

\il<«re  UQ  .Uj  ,U2  ,  are  continuous  twice-differentiable  functions 
of  the  normalised  time  t  and  that  0<  w  <  1  is  the  perturbation 
parameter.  The  initial  conditions  obtained  from  equations  (2)  and 
(3)  in  terms  of  U  are 


UQ(0)  =1,  Ui(0)  =0,  U2  (0)  =0 


NT  - 
■> 
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(4) 
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Substituting  the  expression  for  U  and  equating  coefficients  of  u> 
to  zero  for  r  =  0,  1,  2  one  obtains: 


U  -  i-ajU  =  0 
o  4*o 


Ui  -  ftlVi  =  0 


U2  -  laiU2  «e  2  T 


For  equation  (5)  the  solution  is 


1  1 
-zaix  -5*1* 

U  =  A  e2  +  Be  2 
00  o 


where  A^  and  B  are  constants  to  be  determined  by  the  initial 
00 

coditions,  then 


0*1  T 

Do  =  e* 


Similarly  for  equation  (6)  the  solution  is 


Bl  =  Aje' 


2  +  B,e  2 


where  Ai  and  Bi  are  constants,  and  substituting  the  initial 
conditions  one  obtains 


2III 

Uj  =  - ~r=(  e 

aiII0  JlT 


^  ). 
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Finally  for  equation  (7)  the  solution  is 


1  1  3 

-^ai-t  !  “2alx 

U2  =  A2e  +  B2e  +  (  )e 


where  A2  and  B2  are  constants  to  he  determined  from  the  initial 
conditions,  thus 


”2  =  Tzf  e2 


alT  -paiT  -paxx 

-  2e  2  +  e  2  ) 


(10) 


The  final  solution  for  the  second  invariant  once 
substituting  the  perturbation  parameter  to  by  the  related  value  X 
in  the  expansion  for  U  and  transforming  to  the  original  variable 
II  is: 


2X111.  X2II2  -2aXT 
II  a  (  II  +  — — °  +  )e 


The  final  solution  for  the  third  invariant  is  obtained  by  using 
equation  (1) 


j  XII2  -3aix  XII2  -4a!xj 

111  ■  ino!<  1  +2^rrr  )e  -  izrm-*  i 


(12) 


Authors 


Q 


Tucker  &  Reynolds 

(1968) 


Gence  &  Mathieu 

(1980),  o=0 


le  Penven  et  al, 
(1984)  .(1) 


-0.03213 

0.0175 

0.21 

0.8 

0.7 

0.03502 

0.1157 

0.37 

0.8 

0.7 

0.03418 

0.1191 

0.28 

0.7 

0.7 

0.01729 

0.0775 

0.34 

0.2 

0.5 

0.02487 

0.02251 

0.61 

0.7 

0.4 

-0.4248 

0.03986 

0.5 

0.4 

0.5 

0.  6 


0.4 


0.  2 


0.  0 


ft,® 


\  *  ;v 
t  *  *  %  ' 

i^o 


0.  00 


0.  06 


0.  12 


0.  18 


0.  24 


0.  30 


NORMALISED  TIME 


V-r. 


FIGURE  1 


COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 
COMPONENTS  IN  THE  UBEROI  EXPERIMENT  (1956)  AS  A 
FUNCTION  OF  THE  NORMALISED  TIME.  LINEAR  SOLUTION. 
(CONTRACTION  RATIO  4:1.  R..  =  3710:  DATA: 
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COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 
COMPONENTS  IN  THE  UBEROI  EXPERIMENT  (1956)  AS  A 
FUNCTION  OF  THE  NORMALISED  TIME.  LINEAR  SOLUTION. 
(CONTRACTION  RATIO  *1:1.  R..  =  6150;  DATA: 
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FIGURE  3 


COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 
COMPONENTS  IN  THE  UBEROI  EXPERIMENT  (1956)  AS  A 
FUNCTION  OF  THE  NORMALISED  TIME.  LINEAR  SOLUTION. 
(CONTRACTION  RATIO  4:1,  R„  =  12300;  DATA: 
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FIGURE  4  RELATION  BETWEEN  THE  CONSTANTS  CONNECTED  TO  THE 

LINEAR  AND  NON-LINEAR  PART  OF  THE  PRESSURE-STRAIN 
TERM  FOR  AXISYMMETRI C  FLOWS,  UBEROI  EXPERIMENT 
(1956). 
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FIGURE  5  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 

COMPONENTS  AS  A  FUNCTION  OF  THE  NORMALISED  TIME, 
WITH  VALUES  OBTAINED  FROM  THE  LINEAR  AND  NON-LINEAR 
ANALYTICAL  SOLUTIONS  USING  DIFFERENT  COMBINATIONS 
OF  ci  and  A.  (UBEROI  EXPERIMENT  (1956)  CONTRA CTIO 
RATIO  4:1,  R„  =  3710,  IIKO,  LINEAR - ,  NON¬ 
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FIGURE  6  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 

COMPONINTS  AS  A  FUNCTION  OF  THE  NORMALISED  TIME, 
WITH  VALUES  OBTAINED  FROM  THE  LINEAR  AND  NON-LINEAR 
ANALYTICAL  SOLUTIONS,  USING  DIFFERENT  COMBINATIONS 
OF  Ci  and  A.  (UBEROI  EXPERIMENT  (1956),  CONTRACTION 
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FIGURE  8(a) 
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FIGURE  8(b)  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 
TENSOR  AS  A  FUNCTION  OF  THE  NORMALISED  TIME  IN  THE 
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FIGURE  10(a)  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 

TENSOR  (PRINCIPAL  COMPONENTS)  AGAINST  THE 
NORMALISED  TIME  IN  THE  GENCE  AND  MATHIEU  EXPERIMENT 
(1980),  WITH  PREDICTIONS  OBTAINED  FROM  THE  LINEAR 
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FIGURE  10(b)  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 
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FIGURE  11(a)  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 

TENSOR  (PRINCIPAL  COMPONENTS)  AGAINST  THE 
NORMALISED  TIME  in  THE  GENCE  AND  MATHIEU  EXPERIMENT 
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FIGURE  11(b)  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 

TENSOR  (PRINCIPAL  COMPONENTS)  AGAINST  THE 
NORMALISED  TIME  IN  THE  GENCE  AND  MATBIEU  EXPERIMENT 
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1.1(b)  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 
TENSOR  (PRINCIPAL  COMPONENTS)  AGAINST  THE 
NORMALISED  TIME  IN  THE  GENCE  AND  MATHIEU  EXPERIMENT 
(1980),  WITH  PREDICTIONS  OBTAINED  FROM  THE  LINEAR 
AND  PERTURBATION  SOLUTIONS.  THE  INITIAL  NORMALISED 
TIME  HAS  BEEN  SHIFTED  TO  THE  POINT  WHERE  THE 
ANISOTROPY  STARTS  TO  DECAY  SMOOTHLY. 
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FIGURE  12(a)  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 

COMPONENTS  AGAINST  THE  NORMALISED  TIME  IN  THE  LE 
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FIGURE  12(b)  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 
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FIGURE  15(a)  COMPARISON  OF  THE  EVOLUTION  OF  THE  REYNOLDS  STRESS 
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AND  STIRRED  REACTOR  COMBUSTOR  MODELLING  PROCEDURES 
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SUMMARY 

The  objective  of  this  study  is  to  extend  the  scope  of  finite  difference 
modelling  of  gas  turbine  combustors  so  that  the  maximum  air  loading  (and  also 
the  rich  and  lean  limits)  can  be  predicted.  Since  the  finite  difference 
procedure  is  not  suitable  for  the  direct  computation  of  blow  off  limits,  it 
has  been  adapted  to  calculate  the  residence  time  distribution  from  which  the 
equivalent  stirred  reactor  of  the  primary  stabilisation  zone  can  be  deduced. 
Well  established  stirred  reactor  modelling  can  then  be  used  to  compute  the 
combustion  stability  loop  and  other  features  such  as  pollutant  production 
and  combustion  efficiency.  In  this  paper,  a  pulsed  mercury  tracer,  which  was 
detected  optically,  is  used  to  measure  the  residence  time  distribution  and 
this  result  is  compared  to  that  computed  from  our  finite  difference  code. 

In  the  computation  the  turbulent  trajectories  of  a  large  number  of  neutrally 
buoyant  particles  are  used  to  represent  the  tracer  pulse.  Good  agreement  is 
found  between  the  measured  and  the  computed  response,  and  the  subsequently 
deduced  maximum  air  loading  of  the  combustor  is  consistent  with  the  measured 
value  for  a  variety  of  fuel  nozzles. 

INTRODUCTION 

Gas  turbine  combustors  involve  the  highly  complex  simultaneous  processes 
of  three-dimensional  flow,  droplet  evaporation,  gas  phase  mixing,  radiation 
and  chemical  kinetics.  The  complexity  of  these  interacting  processes  is 
such  that  most  current  gas  turbine  combustor  designs  depend  largely  on 
empirical  correlations.  In  recent  years,  thanks  to  the  rapid  development 
of  computers  and  concurrent  research  into  the  fundamentals  of  the  various 
aspects  of  combustion  and  turbulent  flow,  mathematical  models  of  a  fundamental 
nature  have  started  to  emerge.  Although  not  yet  sufficiently  reliable 
design  tools,  these  models  have  proven  to  be  very  effective  in  analysing 
existing  designs  and  it  is  likely  that  they  will  progressively  replace  the 
cut-and-try  approach  to  combustor  design. 


This  paper  addresses  the  problem  of  closing  the  gap  between  two 
different  approaches  to  mathematically  modelling  the  behaviour  of  practical 
gas  turbines  and  ramjet  combustion  systems.  These  approaches  are:- 

a)  the  finite  difference  method  of  evaluating  the  governing  differential 

equations  ,  and 

b)  the  stirred  reactor  network  method  which  evaluates  the  simultaneous 

governing  equations  for  the  combustor  represented  by  an  interconnected 

2 

set  of  individual  reactors  . 

In  recent  times  the  finite  difference  approach  has  gained  rapidly  in 
popularity  since  it  can  resolve  details  of  the  flow  and  temperature  fields 
which  are  important  to  the  turbine  designer.  Nevertheless  there  are 
significant  features  of  combustor  operation,  such  as  the  prediction  of  the 
rich  and  weak  combustion  stability  limits  as  a  function  of  pressure, 
temperature  and  flow,  which  are  both  difficult  and  expensive  to  compute  via 
the  finite  difference  procedure.  Fortunately,  the  reactor  network  approach 
is  well  suited  to  this  latter  task  and  the  two  approaches  should  be  regarded 
as  complementary  rather  than  as  competing  alternatives. 

For  any  given  combustor  input  pressure  and  temperature,  a  maximum 
combustor  throughput  is  obtained  when  the  rich  and  lean  limits  narrow  until 
they  close  the  stability  loop.  This  maximum  attainable  combustion  intensity 
for  the  whole  chamber  is  an  important  design  parameter  since  it  determines 
the  maximum  operational  altitude  of  the  engine.  Similarly,  since  the 
operation  of  almost  any  kind  of  combustor  at  mixtures  richer  than  stoichio¬ 
metric  results  in  severe  losses  in  efficiency,  we  are  particularly  interested 
in  the  lean  lirjLt  of  combustion  since  this  affects  the  turn  down  range  of  the 
gas  turbine  or  ramjet  engine.  The  difficulty  in  using  finite  difference 
modelling  procedures  for  calculating  the  stability  limits  arises  partly  from 
the  form  of  the  differential  equations.  These  are  stiff  non-linear  equations 
and  their  use  in  searching  for  the  flame  blow-off  loop  conditions  not  only 
involves  long  computation  times  but  is  also  complicated  by  the  need  to 
distinguish  against  numerical  instabilities  when  computing  close  to  a  real 
flame  instability.  The  stirred  reactor  model  is  therefore  potentially  more 
suitable  for  the  accurate  prediction  of  chemical  kinetic  effects  on  combustion 
efficiency  and  pollution,  maximum  throughput  and  stability  limits. 

APPROACH 

To  be  able  to  set  up  the  reactor  network,  knowledge  of  the  various 
reactor  volumes  and  their  sequence  is  necessary.  It  is  the  object  of  this 
study  to  demonstrate  that  one  can  define  these  volumes  by  finite  difference 


modelling  which  involves  the  solution  of  a  large  set  of  coupled  partial 
differential  equations  for  the  transport  of  mass,  momentum,  energy  and  species, 
supplemented  by  appropriate  models  of  turbulence,  turbulence-chemistry 
interactions,  droplet  combustion  and  evaporation  and  radiation. 

The  field  distribution  of  turbulent  dissipation  obtained  by  the  solution 
of  a  partial  differential  transport  equation  can  also  be  utilized  to  define 
the  stirred  and  plug  flow  regions  with  the  degree  of  stirring  or  mixing  rate 
being  determined  from  the  total  dissipation  rate  within  the  reactor. 

To  be  able  to  achieve  these  ends,  experimental  information  on  the 
residence  time  distributions  (RTD)  in  actual  combustors  is  extremely  helpful 
in  not  only  determining  the  reactor  volumes  directly,  but  also  to  validate 
the  more  fundamental  finite-difference  models  if  these  are  to  be  used  with 
confidence  in  the  future.  A  natural  experimental  monitor  for  the  predicted 
residence  time  distribution  is  the  measured  residence  time  distribution 
function  which  is  a  function  of  both  the  input  and  the  output  statistics  and 
which  for  a  selected  input  position  of  tracer,  can  be  expected  to  provide  the 
details  ox  the  history  of  fluid  parcels. 

In  general,  the  RTD  consists  of  a  dead  time  corresponding  to  the  plug 
flow  component  and  an  exponential  response  corresponding  to  the  stirred 
component.  Unfortunately,  if  the  combustor  flow  field  is  modelled  by  a 
network  of  PSR-PFR  elements,  the  overall  RTD  does  not  necessarily  confirm 
the  uniqueness  of  the  model  configuration.  However,  this  problem  can  be 
overcome  by  determining  the  RTD  between  various  locations  within  the  chamber. 

In  this  study,  an  effort  has  been  made  to  measure  the  RTD  functions  in 
a  cold  flow  combustor  using  mercury  vapour  as  a  tracer.  Concurrently  a 
predicted  three-dimensional  velocity  field  and  the  associated  turbulence 
quantities  have  been  utilized  to  simulate  the  experiments  numerically  on  a 
digital  computer.  Detailed  description  of  the  experimental  method  and  the 
theoretical  approach  is  provided  in  the  following  sections. 

EXPERIMENT 

The  geometry  used  for  the  study  was  a  research  Lycoming  gas  turbine 
combustion  can  on  which  much  experience  has  been  gained  over  a  number  of 
years.  This  is  as  shown  in  Fig.  1  and  consists  of  an  upstream  tangential 
air  entry  (7.8%)  around  the  axial  fuel  inlet  and  thereafter  three  groups  of 
six  circumferentially  distributed  holes  at  positions  downstream  representing 
primary  (25.5%),  secondary  (29.9%)  and  dilution  (36.8%)  air  inlet  ports. 


C-; 


r. 


1 


* 


and  fuel  evaporation  patterns.  In  this  procedure ,  individual  droplets  with 
initial  size  and  velocity  are  tracked  through  the  combustion  chamber  as  they 
experience  the  predicted  local  flow  velocities.  At  each  location  a  turbulent 
contribution  is  also  included  whose  magnitude  is  based  on  the  random  choice 
of  a  normally  distributed  variable  simulating  Gaussian  turbulence,  and  whose 
duration  is  based  on  an  eddy  decay  model.  For  a  statistically  large  number 
of  droplets,  this  has  yielded  useful  information  on  evaporation,  droplet 
escape  and  choice  of  fuel  nozzle  for  particular  geometries. 

In  a  similar  manner,  such  a  procedure  may  be  used  for  the  converged 
flowfields  with  a  neutrally  buoyand  fluid  tracer  element  and  if  sufficient 
elements  are  tracked  from  the  fixed  input  position  to  the  exit  plane,  the 
resulting  number  versus  time  is  the  residence  time  distribution  function, 
which  fully  embodies  the  effects  of  broadening  due  to  turbulence  and 
recirculation. 

SIMULATION 

The  kinematic  relations  which  describe  the  motion  of  a  neutrally 
buoyant  particle  in  a  random  turbulent  flow  field  can  be  written  in  a 
cylindrical  polar  system  of  co-ordinates  as 


dZ 

dt 


$  and 


& 


(1) 


where,  t  is  the  time,  u,  v  and  w  are  the  instantaneous  velocity  components 
in  the  Z,  r  and  6  directions  respectively.  The  solution  of  the  above 
stochastic  ordinary  differential  equations  requires  knowledge  of  the 
turbulence  field  in  a  Lagrangian  frame  of  reference.  Unfortunately,  such 
information  is  not  directly  available.  However,  it  is  possible  to  work  out 
the  simulation  using  Eulerian  information.  The  instantaneous  velocities  in 
the  above  equations  can  be  decomposed  into  a  mean  and  a  fluctuating  part  as 

i 

u.  =  u.  +  u. 

Ill 

I 

where,  u^  is  the  time  averaged  velocity  component  and  u^  is  the 
instantaneous  fluctuation.  The  time  averaged  part  of  the  velocity  field 
can  either  be  measured  or  obtained  from  the  solution  of  the  governing 
equations  of  the  conservation  of  mass  and  momentum  supplemented  by  a  suitable 
turbulence  model.  The  instantaneous  fluctuations,  on  the  other  hand,  can 
be  related  to  the  local  kinetic  energy  of  turbulence  K,  under  the  assumption 
of  isotropy  and  Gaussian  probability  density,  as 
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where  +  is  a  normally  distributed  random  variable.  For  isotropic  homogeneous 

turbulence  the  fluctuating  velocities  persist,  on  the  average,  for  a  time 

period  equal  to  the  Lagrangian  integral  time  scale  tt ,  which  is  related  to 

"  6 

the  kinetic  energy  of  the  fluctuating  motion  K  and  its  dissipation  rate  e  by 


=  0.16  - 
E 


(3) 


The  flow  field  is  obtained  in  this  study  from  a  finite  difference 
solution  of  the  equations  of  mass  and  momentum  conservation  in  3-dimensions. 

The  turbulence  is  modelled  by  the  widely  used  k-e  model,  which  entails  the 
solution  of  transport  equations  for  the  kinetic  energy  of  turbulence  and  its 
dissipation  rate.  Full  details  of  these  calculations  can  be  found  in  Ref.  (1). 

The  random  trajectory  of  each  particle  is  calculated  in  the  following 
manner. 

1.  Given  the  location  of  the  particle  at  t®tQ  the  values  of  the 
instantaneous  velocity  fluctuations  and  the  time  during  which  these 
exist  are  obtained  from  eqns  (2)  and  (3).'  It  is  assumed  that  the 
mean  velocity  components  remain  constant  within  each  finite  difference 
cell  that  the  particle  is  traversing  and  change  abruptly  from  cell  to 
cell.  No  effort  is  made  to  obtain  a  local  value  by  interpolation  as 
this  would  be  inconsistent  with  the  approximate  nature  of  the  simulation. 

2.  The  kinematic  relations  (1)  are  solved  for  a  time  step  which  is  the 
smaller  of  x^  and  the  time  required  for  the  particle  to  cross  the 
current  finite  difference  cell. 

3.  At  each  cell  crossing  the  appropriate  values  of  the  mean  velocities 

are  inserted  into  eqns  (1).  New  instantaneous  fluctuations  are 

obtained  from  eqn  (2)  and  x^  from  Eq.  3  when  the  particle  has  travelled 

for  a  time  interval  equal  to  xT. 

L 

4.  Steps  (1)  to  (3)  are  repeated  until  the  particle  leaves  the  combustion 
chamber. 

The  simulation  technique  is  equivalent  to  solving  the  time  dependent 
transport  equation  of  a  passive  scalar,  which  would  require  considerable 
computing  effort  especially  in  3-D  environments. 
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RESULTS  AND  DISCUSSION 


I 


Using  the  experimental  technique  described  above  a  set  of  measurements 

of  the  RTD  were  recorded.  Individual  absorption  traces  exhibit  variation 

due  to  background  noise  and  variation  of  the  mixing  history  for  each  pulse 

of  tracer.  Although  these  can  be  very  useful  for  obtaining  gross  mixing  proper- 

3 

ties  or  for  extracting  empirical  constants  for  mixing  models  their  average 
effects  have  particular  significance  when  comparing  with  results  from  funda¬ 
mental  flow  and  turbulence  modelling.  Each  experiment  thus  consisted  of 
an  average  of  SO  such  individual  pulse  measurements  and  as  a  guard  against 
absolute  pulse  to  pulse  variation,  each  signal  was  normalised  to  a 
predetermined  area. 

The  resulting  trace  approximated  to  the  smoothed  pulse  response  of 
the  system  for  the  particular  input  position  considered  and  an  example  is 
given  in  Fig.  3. 

Five  trace  input  positions  were  investigated  both  experimentally  and 
computationally.  All  tracer  output  measurements  were  made  at  the  exit  from 
the  chamber  and  moving  progressively  upstream,  the  input  positions  were:- 

a)  through  the  secondary  port, 

b)  on  the  axis  between  the  primary  and  secondary  ports, 

c)  on  the  axis  at  the  primary  port  position, 

d)  through  the  primary  port, 

e)  on  the  axis  at  the  fuel  nozzle  position. 

For  the  computation  results,  about  6000  'particles'  were  tracked  for 
each  inlet  location  to  obtain  an  adequate  statistical  sample. 

The  residence  time  distribution  functions  both  measured  and  simulated 
are  shown  in  Fig. 4, 5.  The  agreement  between  the  experimental  and  simulated 
residence  time  distributions  are  remarkably  good  for  the  three  shorter 
distances,  however  the  discrepancy  increases  for  the  two  longer  distances. 

Errors  are  likely  to  be  due  to 

i)  the  use  of  a  computed  flow  field  instead  of  a  measured  one, 

ii)  to  the  use  of  the  k-e  model  results  to  approximate  the  Lagrangian 
information, 

iii)  to  neglect  of  the  finite  duration  of  the  input  pulse,  and 

iv)  the  possible  existence  of  some  computational  dead  time  on  the  small 
axial  finite  difference  cell. 
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Nevertheless,  the  agreement  is  generally  satisfactory  and  confirms  the 
validity  of  this  method  of  extracting  residence  time  distribution  information 
from  the  finite  difference  modelling. 

The  next  step  is  to  convert  the  tracer  information  to  dead  time  (plug 
flow)  and  stirred  reactor  time  constant  components.  In  particular  we  wish 
to  characterize  the  primary  zone  since  this  controls  the  combustion 
stability.  As  pointed  out  above,  the  stirred  reactor  regions  can  be 
identified  readily  from  the  computed  regions  of  recirculating  flow  and  high 
turbulence  intensity  or  dissipation  . 

In  the  cases  of  tracer  injection  through  the  secondary  port,  or  on  the 
axis  between  the  primary  and  secondary  ports,  inspection  of  the  responses  of 
Figure  4  shows  very  similar  residence  time  functions.  These  locations 
represent  the  exit  from  the  primary  zone,  thereiore  to  characterize  the 
primary  zone  from  the  overall  response  function,  we  must  first  analyse  the 
response  function  from  either  of  these  locations. 

The  integrated  output  signals  for  tracer  injection  through  the  secondary 
zone  and  primary  zone  are  plotted  on  semi-log  co-ordinates  in  Fig.  6.  From 
this  figure  it  is  clear  that  the  response  for  the  injection  into  either  zone 
consists  of  a  dead  time  followed  by  an  exponential  decay.  The  response  of 
the  primary  zone  alone  can  be  obtained  by  difference.  The  values  are  as 
follows  :- 

Dead  time  Stirred  reactor 

(ms)  time  const. (ms) 

Injection  into  secondary  zone  5.1  9.8 

Injection  into  primary  zone  9.3  16.8 

Primary  zone  alone  4,2  7.0 

The  primary  zone  transfer  function  io  thus: 

G(s)  =  exp(-0. 0042s) /(0.007s+l) 

The  effective  volume  of  the  primary  zone  PSR  can  be  deduced  directly  from 
this  time  constant  since  the  volumetric  flow  rate  into  the  primary  zone  is 
known  to  be  0.03m  /s.  Hence  the  effective  primary  zone  volume  is  0.00021  mJ. 
This  result  is  consistent  with  the  value  estimated  by  inspection  of  the 
volume  of  the  chamber  and  the  computed  size  of  the  recirculation  zone.  It 
therefore  confirms  that  the  whole  of  this  zone  behaves  as  a  single  stirred 
reactor  and  this  supports  the  use  of  the  PSR  model  for  predicting  combustor 
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behaviour.  In  order  to  calculate  the  lean  limit  of  the  combustor  at  any 
particular  input  pressure  and  temperature ,  and  fuel  spray  size  distribution, 
we  would  have  to  know  the  local  equivalence  ratio.  Fortunately,  this 
information  is  available  from  the  finite  difference  model,  and  due  account 
can  be  taken  of  the  important  regions  where  the  mixture  is  locally  stoichio¬ 
metric.  Hence  the  two  modelling  procedures  have  complementary  utility  to 
the  designer. 

The  maximum  air  loading  ip  of  the  combustor  is  attained  when  the 

equivalence  ratio  of  the  whole  primary  zone  is  unity.  The  actual  overall 

fuel/air  ratio  at  which  this  is  obtained  depends  on  the  droplet  size 

distribution  and  the  operating  conditions,  however,  the  assumption  of  a 

stoichiometric  PSR  allows  the  well  established  relations  for  stirred  reactors 

7  8  9 

to  be  used.  In  the  case  of  kerosene  fuel,  the  appropriate  relation  is  *  * 


_  &  _  K  exp(-E/RT  )<l-n  )°‘75(l-*n  ) 

r  “*  “ L  T  jc“  ^  C  C 

VP  T_1.25  0.25 

R  *  nc 

where  the  equivalence  ratio  $  should  be  taken  equal  to  unity,  and  in  SI 

units,  E/R  =  21,000  and  K  *  8400.  V  =  volume,  P  =  pressure,  T  =  reactor 

temperature ,  n  *  combustion  efficiency  and  m  =  air  mass  flow.  Thus  for 
C  5  2 

the  cease  in  question,  at  10  N/m  ,  tj>  «*  1,  and  Tq  =  350  K,  Ref.  7  Fig  5.7 
shows  that  fic  “  0.78  at  the  point  of  maximum  reaction  rate  in  a  stirred 
reactor.  Ref.  7  Fig.  2.20  gives  an  adiabatic  temperature  rise  of  1850  K  at 
these  conditions ,  hence  the  primary  reactor  temperature  is  350+0. 78xl850z1800. 
The  maximum  air  loading  is  therefore  predicted  to  be: 


=  fc  =  8400  exp(-21, 000/1800)  Q.»22 
max  vP1*75  18001’25  0.78 

=  5.6xl0-7 

.*.m  =  5.6xl0~7  x  0.00021  x  (105)1* 75 

max 

=  0.066  kg/s. 

Thus  the  total  chamber  flow  at  blow  off  0.2  kg/s.  The  validity  of  this 
result  is  illustrated  in  Figs  7a)  and  b)  in  which  measured  stability  loops 
for  this  combustor  are  plotted  for  a  variety  of  fuel  nozzles  (Details  of  the 
experiment  are  given  in  Ref. 10).  It  is  clear  that  the  maximum  air  loading 
for  all  the  fuelin'  systems  tested  occur  close  to  the  expected  value. 

It  is  also  apparent  that  the  technique  presented  here  can  be  readily 
extended  to  the  prediction  of  rich  and  lean  combustion  limits,  however  this 
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Figure  1.  Gas  turbine  can  geometry  and  dimensions. 

Figure  2.  Example  of  individual  tracer  pulse  response. 

Figure  3.  Tracer  response  signal  obtained  by  averaging  50  individual 

single  shot  responses. 

Figure  4.  Comparison  of  observed  tracer  response  signal  and  calculated 

residence  time  functions  for  tracer  input  to  the  primary 
zone  of  gas  turbine  can.  - -  =  observed  response , 

•  =  calculated  values. 

Figure  5.  Comparison  of  observed  tracer  response  signal  and  calculated 

residence  time  functions  for  tracer  input  to  the  secondary 
zone  of  gas  turbine  can.  * — - -  observed  response , 

•  =  calculated  values. 

Figure  6.  Log  of  integrated  pulse  response  (step  response)  vs  time  for 

combined  primary  zone  responses  and  combined  secondary  zone 
responses.  O  =  secondary  zone  response,  A  *  primary  zone 
response. 

Figure  7.  a)  Stability  loop  for  pressure  jet  atomisers  -  kerosene  fuelled 

b)  Stability  loop  for  air  blast  atomisers  -  kerosene  fuelled 
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RESIDENCE  TIME  FUNCTION 

FOR  AXIAL  TRACER  INPUT  AT  PRIMARY  PORT  POSITION 


FUEL  RATIO 


A  NO  PHOT  FLOW 
o  PILOT  JET.  791  kPo 
o  PILOT  JET.  543  kPa 


FIGURE  7 


Introduction 


The  measurement  of  the  size  distribution  of  droplets 
and  particles  is  a  common  problem  in  the  current  technology. 
In  many  processes  measurement  and  control  of  particle  size  is 
of  great  importance.  An  optical  technique  initially  developed 
in  this  department  and  later  in  conjunction  with  Malvern 
Instruments  Ltd,  makes  on-line  particle  sizing  possible.  The 
technique  uses  the  diffraction  pattern  formed  when  particles 
are  illuminated  by  a  parallel  beam  of  monochromatic,  coherent 
light  to  determine  the  particle  size  distribution. 

One  6f  the  major  limitions  of  this  method  is  that  for 
reliable  results  not  more  than  fifty  percent  of  the  incident 
light  should  be  scattered  by  the  particle  field.  It  has  been 
speculated  that  above  this  limit  the  effect  of  multiple 
scattering  of  light  by  the  particles  will  have  a  great 
influence  on  the  results. 

However  in  many  problems  such  as  those  encountered  in 
particle  sizing  in  diesel  injectors  and  heavy  fuel  oil 
atomizers  this  limit  is  usually  exceeded.  A  technique  has  been 
developed  which  predicts  the  effect  of  multiple  scattering  on 
the  diffration  pattern.  For  distributions  which  can  be 
described  by  simple  functions  such  as  Rosin-Rammler  and 
Log-Normal  a  set  of  correction  equations  have  been  derived  to 
compensate  for  the  effects  of  multiple  scattering.  As  an 
illustration  of  this  approach  results  are  presented  for  a 
commercial  diesel  injector  which  gives  up  to  90%  obscuration. 
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When  a  spherical  particle  is  illuminated  by  a  parallel 
beam  of  monochromatic,  coherent  light  a  Fraunhofer  diffraction 
pattern  is  formed,  superimposed  on  the  geometrical  image,  this 
pattern  is  large  compared  with  the  image.  If  the  particle 
field  is  placed  within  the  focal  length  of  a  suitable  Fourier 
transform  lens  then  the  undiffracted  light  is  focussed  to  a 
point  on  the  axis  of  the  lens  and  the  diffracted  light  forms  a 
far  field  Fraunhofer  pattern  of  rings  around  the  central 
spot (Fig. 1) .  Movement  of  the  particles  does  not  cause  movement 
of  the  diffraction  pattern  since  light  diffracted  at  an  angle, 
6,  will  give  the  same  radial  displacement  in  the  focal  plane 
irrespective  of  the  particle  position  in  the  illuminating 


beam. 


The  theory  and  analysis  of  light  scattering  by  single 


spherical  particle  can  be  found  in  standard  texts  on  optics1 
and  has  been  briefly  discussed  in  a  paper  by  Swithenbank  et 
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The  diffraction  pattern  for  a  single  sphere  (or  a  field 
of  monosize  spherical  particles)  is  shown  in  figure(2).  The 
light  intensity  is  maximum  at  the  centre  of  the  diffraction 
pattern  and  oscillates  with  strongly  decreasing  amplitude  as 
the  radius  increases. 

The  diameter  of  the  diffraction  pattern  is  inversely 
proportional  to  the  particle  diameter  and  the  smallest 
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particle  that  can  be  measured  is  of  the  order  of  one  micron 
since  the  particle  size  must  be  greater  than  the  wavelength  of 
the  illuminating  laser  light (wavelength=0 .6328um) . 

As  the  diffraction  pattern  is  unique  for  a  given 
particle  size  or  a  given  mixture  of  particles  then  the  light 
energy  pattern  can  be  used  to  deduce  the  size (or  size 
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distribution)  of  the  particle(s). 

For  a  monosize  distribution  of  particles  of  size  a,  the 
light  energy  contained  within  a  circle  of  radius  s  in  the 
focal  plane  of  a  lens  is  given  by:- 


L  =  1  -  J 


2  2iras 


2  2iras 
1  Xf 


Thus  the  light  energy  within  eny  ring  in  the  focal 
plane, bounded  by  radii  s,and  s2is  given  by:- 

Ls  c  =  e  jjT0  (Kas2/f)  +  (Kas^/ f )  -  J^CKas^/f)  -  J^Kas^f)  .  (2) 

Where  £  is  the  energy  falling  on  the  particle  which  is 
proportional  to  the  cross-sectional  area  of  the  particle, 
therefore:- 


E  =  C  *  Tra  .N 


And  N  is  given: 


Thus: - 


N  =  3w/ .  3 

4ira  p 


E  =  C"  - 


....  (4) 


Where  W  is  the  weight  fraction  of  particles  of  radius  a.  C' 
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and  C"  are  constants  which  depend  on  the  power  of  the  laser 
used,  p  is  the  density  of  the  particle.  Figures  (3)and(4)  are 
plots  of  equations  (1)  and  (2) .  If  instead  of  a  single 
particle  a  collection  of  particles  of  different  size  is 
considered  then  the  light  energy  falling  on  any  ring  in  the 
focal  plane  is  the  sum  of  the  contribution  from  individual 
particles: - 
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w. 
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j  =  C0  \ 

sl’s2  2  /  77 

L—l  1 

i=l 


0  Ka.  .  Ka. 

T2/  1  V  .  j2  /  1  \ 

Jo{~  si}  +  Ji  (1“  si) 


0  Ka.  _  Ka. 

Jo  (_r  S2)-Ji <-r  V 
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(6) 


Where: 

K  =  2v/\ 


(7) 


If  a  particle  size  distribution  is  classified  into  a  number  of 
size  ranges  then  it  is  preferable  to  make  measurments  at  radii 
where  the  first  maximum  of  the  diffracted  energy  distribution 
occurs.  This  location  is  given  by:- 

2iTas/xf  =  1.375  . (8) 

If  a  derector  is  \ised  which  is  divided  into  a  set  of  circular 
rings  then  each  of  these  rings  will  define  a  characteristic 
particle  size  range.  The  total  light  energy  distribution  is 
the  sum  of  the  product  of  the  energy  distribution  for  each 
size  range  and  the  weight  fraction  in  that  range.  This  can  be 
expressed  in  the  form  of  a  matrix  equation:- 

L(I)  =  W(J)  T(I , J)  . (9) 
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Where  L(I) 


and  W ( J )  are  the  light  energy  and  weight 
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distribution  respectively  and  T(I,J)  contains  the  coefficients 


which  define  the  light  energy  distribution  for  each  particle 


size  range.  The  weight  distribution  is  deduced  in  such  a  way 


that:  - 


I(L(J)  -  W(J)  T  (I,J))2  is  a  minimum 


The  ratio  of  the  light  intensity  measured  at  the  centre 


diode,  before  and  after  the  particle  field  is  placed  in  the 


laser  beam,  gives  the  fraction  of  the  light  transmitted  by  the 


particles.  The  transmission  is  related  to  the  total  projected 


area  of  the  particles  by  the  Beer-Lambert  law: 


=  -  tS, 


Where  I  and  I0are  the  light  intensities  with  and  without  the 


sample  respectively,  £  is  the  optical  path  length  and  t  is 


given  by:- 


t  =  2  \  N.A. 

/  .  1  1 
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In  deducing  the  size  distribution  from  the  light  energy 


distribution  several  distribution  functions  have  been  used. 


Normal,  Log-Normal,  Rosin-Rammler  and  Model  Independent  are 


some  examples.  The  distributions  employed  in  this  study  are 


mainly  Rosin-Rammler  and  Log-normal  which  are  both  two 


parameter  forms. 


The  Rosin-Rammler  distribution  which  is  of  particular 


intarest  in  spray  studies  is  defined  by*:- 


Y  =  100(1“  Exp  (~x/x)N) 
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Where  Y  is  the  percentage  of  particles  with  diameter  less  than 
X,  X  is  the  Rosin-Rammler  mean  diameter  and  N  is  the 
Roain-Rammler  exponent , which  is  a  measure  of  the  spread  of  the 
particles  about  the  mean  size  T. 


The  log-Normal  distribution  is  defined  by*:- 


L  =  1 . - 

d  /2  £n(N)  2  £ n  (N)2 


.  (13) 


Where  is  the  weight  frequency  at  size  d,  )C  is  the  geometric 
mean  of  the  distribution  and  N  the  geometric  standard 
deviation. 

The  calculated  size  distribution,  whichever  model  is 
used,  is  reliable  only  if  the  light  obscuration  is  in  the 
range  5-50%.  The  lower  limit  is  the  limit  of  measuring  a 
representative  sample  of  particles  and  the  upper  limit  is  the 
onset  of  significant  of  multiple  scattering.  Within  this  range 
the  sample  measured  is  taken  to  be  representative  of  the 
particle  field  and  the  light  scattered  by  particles  is  not 
further  scattered  by  others. 

Since  the  angle  of  the  scattered  light  varies 
proportionally  with  the  inverse  of  the  particle  diameter,  the 
effect  of  multiple  scattering  is  to  represent  a  bigger 
proportion  of  smaller  particles  than  the  particle  field 
actually  contains. 

The  model  developed  in  this  study  predicts  the  effect 
of  multiple  scattering  on  the  parameters  of  the  distributions 


defined  by  equations  (13)  and  (14). 

In  the  development  of  the  model  the  particle  field  is 


divided  into  slices.  The  following  assumptions  were  made: 

A)  The  light  path  is  divided  into  a  series  of  slices  of 
equal  obscuration  each  slice  scatters  10%  of  the  incident 
light. 

B)  Only  half  of  this  scattered  light  is  forward  scattered 
light (the  scattered  light  as  seen  by  the  detector 
plane-Babinet' s  principle). 

C)  There  is  no  multiple  scattering  within  a  slice. 

D)  The  scattered  light  pattern  leaving  one  slice  is  the 
incident  light  for  the  next. 

The  model  enables  the  determination  of  the  overall 
scattered  light  energy  distribution  for  each  slice.  As  the 
number  of  slices  is  increased  the  fraction  of  the  primary 
laser  beam  that  has  been  scattered  increases  and  hence  higher 
obscurations  can  be  simulated. 

Determination  of  the  light  energy  d.i  stribution  for  the 
first  slice  is  straight  forward  since  the  only  illuminating 
light  is  the  primary  undiffracted  beam.  Equation (6)  is  used 
directly  to  calculate  this  data.  For  the  second  and  all  other 
slices  there  is  the  full  spectrum  of  the  scattered  light  at 
all  angles  as  well  as  the  undiffracted  light  to  be  considered. 

A  numerical  integration  method  is  needed  to  calculate 
these  data.  The  integration  procedure  requires  the  division  of 
the  rings  into  segments.  By  considering  each  segment  as  a 
'point  scatterer'  which  scatters  the  light  that  would  have 
been  incident  on  this  segment  according  to  equation (6)  and 
assumptions  (A)  and  (B) ,  thus  determining  the  light  scattered 
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to  all  other  segments.  Figure (5)  shows  a  typical  integration 
step.  Ring  I  is  the  'scattering  ring'  and  ring  J  is  the 
'receiving  ring'.  The  segment  on  ring  I  is  the  point 
scatterer,  all  the  light  within  this  segment  is  represented  by 
point  S.  The  receiving  segment  R  is  treated  as  an  area. 

Distances  OS,  OA  and  OB  are  known,  the  angle  a  i:~ 
incremented  from  0  to  it  depending  on  how  the  rings  ar  ,* 
divided  into  segments.  Distances  SB  and  SA  can  be  calculated 
using  the  Cosine  equation: - 


(SA) *  = ( O A ) 2+(OS) 2  - 2 ( 0 A )  (OS) Cosa 


(14a) 
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(SB) 2  = (OB) 2+(OS) 2 -2 (OB)  (OS) Cosa 


(14b) 
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In  order  to  simplify  the  geometry  the  point  S  is  transferred 
to  S'  a  point  along  the  axis  OAB  so  that:- 


SP.=S  '  R=  (SA+SB)  /2 


i*  «  * 


Now  using  S'A  and  S'B  in  equation(6)  the  scattered  light 
energy  due  to  S  within  these  radii  can  be  calculated.  Knowing 
the  area  of  the  receiving  segment  and  area  of  the  ring  bounded 
by  S'A  and  S'B  the  fraction  of  the  scattered  light  from  S  to  R 


-  !«*.*'*  •< 
•  V*  *.  -  t. 


is  determined.  This  procedure  is  carried  out  taking  two 
segments  at  a  time  until  the  detector  is  fully  covered. 
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The  total  light  energy  on  any  ring  after  each  slice  is 
the  sum  of  all  the  contributions  made  from  other 
rings (secondary  scattered  components)  plus  the  scattered  light 
from  the  primary  beam (primary  scattered  component)  that  falls 
on  that  ring. 

There  are  a  number  of  points  that  need  to  be  mentioned 
at  this  stage. 

1)  The  rings  on  tne  detector  are  not  of  equal  width. 
The  width  of  the  31  rings  increases  from  0.094mm  to  1.565mm. 
This  led  to  problems  with  the  integration  as  will  be  shown 
later. 

2)  As  a  first  attempt  the  rings  were  divided  into  equal 
area  segments  so  that  the  first  ring  and  the  last  ring  had  4 
and  6972  segments  respectively.  At  this  stage  the  detector  was 
assumed  continuous  and  the  gaps  between  the  rings  were  not 
considered.  This  arrangement  was  found  to  be  totally 
impractical  due  to  the  very  long  computer  time  needed  and  was 
abandoned  in  favour  of  dividing  the  rings  into  equal 
angles (fig. 5) .  This  led  to  a  vast  improvement  in  computer 
time. 

At  this  stage  a  new  method  for  determining  the  amount 
of  scattered  light  within  a  ring  was  developed  which  made  the 
use  of  equation (6)  and  hence  Bessel  functions  redundant  and 
this  also  greatly  reduced  computer  time.  The  method  required 
the  primary  scattered  light  energy  distribution  for  a  given 
sample  of  particles  for  twice  the  diode  dimension.  This  was 
found  either  using  equation (6)  or  by  using  Malvern  software, 
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th  e  latter  option  being  preferred.  A  cubic  spline  fit  was  made 
to  this  data  and  interpolation  was  used  to  find  the  energy 
between  any  two  radii  S' A  and  S'B.  All  the  rest  of  the 
calculation  procedure  remaied  the  same. 

The  reason  that  this  technique  saves  computer  time  is 
due  to  the  large  number  of  times  which  the  Bessel  function  is 
called  and  that  this  function  is  calulated  via  a  ten  term 
polynomial1*  as  opposed  to  the  four  terms  of  the  cubic  spline. 
In  addition  the  Bessel  function  is  called  four  times  per  diode 
ring  and  fifteen  times  per  size  class  whereas  the  cubic  spline 
is  only  called  twice  per  diode  ring. 

As  mentioned  earlier  the  rings  on  the  detector  are  not 
of  equal  width.  Figure (6)  shows  the  detector.  It  was  found 
that  due  to  the  comparatively  ’big'  and  'long'  segments  in  the 
outer  rings  the  energy  on  the  inner  ring  was  grossly 
overestimated  resulting  in  unrealistic  energy  distribution 
curves (fig. 7) . 

Realising  the  cause  of  the  problem  it  was  found  that  it 
was  drastically  reduced  by  dividing  the  bigger  rings  into 
smaller  sub-rings.  This  procedure  was  followed  for  increasing 
numbers  of  subrings  until  no  further  appreciable  change  in  the 
energy  distribution  curve  was  observed.  The  final  arrangement 
consisted  of  having  109  sub-rings  of  comparable  width  so  that 


the  innermost  ring  had  one  ring (sub-ring)  and  the  outer  ring 
was  divided  into  ten  sub-rings.  Table (1)  gives  the  dimensions 
of  the  detector.  This  arrangement  of  equal  angle  segments  and 
109  sub-rings  is  the  one  used  to  derive  all  the  results 


mentioned  in  the  next  section. 

The  energies  calculated  for  all  109  sub-rings  at  the 
end  of  each  slice  were  reduced  to  31  according  as  how  each 
ring  had  been  divided  into  sub-rings.  The  gaps  between  the 
detector  rings  are  all  of  width  0.036mm  which  is  comparable 
with  the  inner  ring  widths  but  negligible  compared  with  the 
outer  ring  widths.  Hence  the  effect  of  the  gaps  was  taken  into 
consideration  which  further  reduced  the  energy  on  the  inner 
rings  in  comparison  to  the  outer  ones. 

Figure (8)  shows  the  final  theoretical  light 
distribution  for  a  given  size  distribution  at  several 
obscurations.  The  peak  of  the  curve  moves  to  the  right  ie  the 
outer  rings  which  represent  smaller  size  classes.  These  curves 
clearly  indicate  the  effect  of  multiple  scattering. 
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Results: 


Theoretical  and  experimental  results  for  a  number  of 
different  fine  particulate  samples  are  explained  and  shown  in 
the  first  part  of  this  section.  Having  justified  the  accuracy 
and  applicability  of  the  model  a  number  of  different 
distributions  are  treated  theoretically  and  correlation 
equations  are  derived. 


A)  four  reference  samples  were  available  for  experiments: 


DVN  Iatex-El9 
DVs  latex-G7 
NBS  SRM  1003 
NBS  SRM  1004 


monosize  13.8  u 
monosize  40.8  u 
range  5.0-30.0  u 
range  37.-105.  u 


The  apparatus  is  shown  schematically  in  figure  (9)  and 
the  cell  used  for  suspension  of  the  samples  is  shown  in  figure 
(10) .  A  very  dense  suspension  of  the  sample  was  made  and  the 
stirrer  speed  was  controlled  and  adjusted  so  that  all  the 
particles  were  in  suspension,  while  not  so  fast  as  to  cause 
cavitation  and  hence  formation  of  bubbles.  A  reading  was  taken 
and  by  removing  one  ml  of  the  suspension  and  replacing  it  by 
one  ml  of  fresh  water  the  concentration  was  reduced,  another 
reading  taken  and  this  procedure  repeated  until  the 
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obscuration  reading  was  less  than  40%.  The  water  used  to  make 
up  the  suspension  had  been  previously  filtered  to  remove 
particles  greater  than  approximately  one  micron.  A  drop  of 
Ncnidet-P42  was  added  as  a  dispersant.  Several  runs  for  each 
of  the  four  samples  were  carried  out  to  ensure  that 
representative  sample  had  been  considered. 

The  model  independent  option  was  used  to  process  the 
data  acquired  for  the  monosize  samples.  The  theoretical  model 
was  used  to  generate  the  light  energy  distribution  for  one  of 
the  monosize  samples (Iatex-El9) .  The  theoretical  and 
experimental  results  are  plotted  in  figures  (11)  and  (12) .  The 
two  plots  show  the  same  trend  as  the  obscuration  is  increased. 
This  clearly  shows  the  applicability  of  the  model  to  monosize 
as  well  as  polysize  particle  fields. 

The  Rosin-Rammler  distribution  was  used  to  process  the 
light  energy  data  and  deduce  X  and  N  for  the  two  SRM  samples. 
Two  correction  parameters  are  defined  as  follows: 


fc-  i*. 


t-JU. 
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C^= (Actual  X)/ (Apparent  X) 


(Actual  N)/ (Apparent  N) 


(16a) 


(16b) 


As  the  obscuration  is  decreased  X  and  N  increase  until,  at 
about  60%  obscuration,  these  values  remain  unchanged  for  any 
further  reduction  in  obscuration. 

Figures  (13)  and  (14)  show  plots  of  and  Cn  .vs. 


I'-TI 


obscuration  for  one  of  the  samples  used(SRM  1003).  The  values 
of  ~X  and  N  at  50%  obscuration  were  taken  to  be  the  true 
representation  of  the  sample  used.  The  model  was  used  to 
generate  the  light  energy  distribution  for  fifty  slices  which 
corresponds  to  approximately  99.5%  obscuration.  These 
theoretical  values  are  also  plotted  in  figures  (13)  and  (14) . 
As  it  can  be  seen  there  is  very  good  agreement  between  the  two 
results. 

Similar  results  were  obtained  for  the  other  sample 
(SRM-1004)  but  for  a  given  obscuration  the  correction  factors 
Cx  and  Cn  are  different  from  the  results  for  the  first  sample 
(SPM-1003) .  This  indicated  that  the  correction  factors  are  not 
only  a  function  of  obscuration  but  also  a  function  of  X  and/or 
N  for  the  sample  considered.  This  point  was  further 
investigated  and  is  discussed  in  the  second  part  of  this 
section. 
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B)  As  mentioned,  the  results  indicate  that  the  correction 
factors  for  X  and  N  are  functions  of  actual  X  and  N  as  well  as 
the  obscuration.  Five  size  distributions  of  approximately 
equal  N(2.3)  and  different  )f(27.-85.  micron)  were  considered. 
For  each  of  these  the  light  energy  distributions  for  up  to 
fifty  slices(99.5%  obscuration)  were  generated  and  using  the 
Malvern  software  these  energy  distributions  were  used  to 
generate  apparent  size  distributions.  For  each  size 


distribution  the  factors  C  and  C  were  calculated  and  plotted 
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a  given  obscuration  show  no  consistant  trend  for  different 
samples  and  that  all  the  correction  factors  for  a  given 
obscuration  are  very  close  to  each  other,  it  was  concluded 

that  C  and  C  are  independent  of  the  actual,  and  hence 

x  n 

apparent,  >f.  The  slight  variation  is  within  the  range  of 
convergence  errors  when  calculating  the  size  distribution  from 
light  energy  distribution. 

Six  size  distributions  of  approximately  equal  1((35. 
micron)  and  different  N{1.4-3.8),  the  range  of  interest  for 
sprays,  were  considered  in  the  same  way  as  the  previous 
investigation,  the  results  are  plotted  in  figures  (17)  and 
(18) .  In  this  case  a  clear  and  consistant  dependence  of  Cx  and 
Cn  for  a  given  obscuration  on  actual  or  apparent  N  is 
illustrated. 

It  was  decided  to  deduce  the  correction  equations  in 
terms  of  obscuration  and  the  apparent  N,  since  this  is  the 
information  that  is  known  when  a  sample  is  measured  and  the 
resulting  light  energy  distribution  processed.  Figures  (19) 
and  (20)  show  these  results  plotted  as  or  Cn  against  the 
apparent  N  for  four  different  obscurations.  Using  multiple 
linear  regression  the  following  correction  equations  were 
found  to  fit  the  theoretical  data  well.  These  equations  apply 
to  Rosin-Rammler  distribution  only:- 


C 

x 


=  1.0  +  (0.036+ .4947 (OB) 
=  1.0  +  (0.035+. 1099(0B) 


8.997.  n(1.9-3.437(OB))+2  ^ 
app 

.  (17) 

8.6501  <0.35-1. 45<0B» 

app 

.  (18) 

These  equations  apply  for  obscurations  in  the  range  65-98%  and 
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for  apparent  N  in  the  range  1.2-3. 8  with  the  given  uncertainty 
limits.  For  obscurations  less  than  65%  the  correction  factor 
is  comparable  with  experimntal  and  convergance  uncertainty. 

The  same  analysis  was  carried  out  for  the  Log-Normal 
distribution.  The  results  are  shown  in  figures  (21)  and  (22) , 
again  the  representative  size  T  decrease  with  increasing 
obscuration (C  increasing).  On  the  other  hand  the  correction 
factor  for  N  decreases,  since  as  the  obscuration  increases  the 
distribution  becomes  broader  and  the  geometrical  standard 
deviation  increases. 

Figures  (23)  and  (24)  are  plots  of  C  and  C  against 
the  apparent  N  for  different  obscuration.  The  Log-Normal 
results  were  not  as  consistent  and  smooth  as  those  for 
Rosin-Rammler ,  this  is  possibly  due  to  the  definition  of  the 
•distribution  and  its  applicability  to  the  problem  since  the 
Log-Normal  distribution  shows  a  long  tail  at  large  sizes. 

Again  multiple  linear  regression  was  used  to  generate 


equations  for  and  Cn  in  terms  of  apparent 


obscuration. 

.  ,,  ,  0,  ,  -4  7. 6371 (OB).  - 

=(1.0+6.87x10  e  )  e 


11  3131 

0. 1395+1. 7843(OB)) 

N 

app 


N  and 

.  (19) 


CL  =(0.96+7.3965xl0_4  e5,6798(OB))-(0. 0129+0. 1415(OB)7*1516)N  .  (20) 

n  app 


These  equations  have  uncertainty  limits  of  ±15%,  for  better 
results  the  following  five  equations  for  and  Cn  can  be 

used.  These  equations  apply  to  obscuration  in  the  range  65-98% 
and  apparent  N  values  of  1. 6-3.0.  The  general  form  of  the 
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equations  are  given  below  and  the  particular  coefficients  are 
listed  in  Table (2). 


C  =A(Exp(-B/N  ) 
x  ^  '  app' 


(21) 


C 

n 


=C° 


D(NaPP> 


(22) 


For  intermediate  obscurations,  linear  interpolation  should  be 
used. 

As  this  work  was  carried  out  with  reference  to  diesel 
sprays  an  example  for  illustration  purposes  is  given  here. 
Figures  (25)  and  (26)  show  the  actual  and  apparent  results  for 
a  diesel  injector.  The  operating  and  experimental  conditions 
are  listed  in  Teble(3).  Figure(27)  is  the  time-obscuration 
plot  indicating  the  duration  of  the  spray. 


9 


Conclusion: 


A  mathematical  model  has  been  developed  to  predict  the 
effect  of  multiple  forward  scattering  in  a  dense  particle 
field  on  the  measured  Rosin-Rammler  and  Log-Normal 
distribution  parameters. 

The  theoretical  results  from  the  model  were  tested 
against  experimental  results  for  equivalent  size 
distributions.  Good  agreement  was  observed  as  a  justification 
of  the  accuracy  of  the  model.  Two  correction  factors  have  been 
defined  and  the  dependency  of  these  factors  on  obscuration  and 
distribution,  parameters  (IT  and  N)  have  been  examined. 
Consequently  a  number  of  correction  equations  have  been 
derived  which  allow  the  applicability  of  this  method  of 
particle  sizing  be  extended  twenty  five  fold  from  an 
obscuration  of  50%  to  that  of  98%.  The  uncertainty  limits  of 
using  these  correction  factors,  especially  for  the 
Rosin-Rammler  distribution,  is  comparable  with  experimental 
and  solution  convergence  uncertainties. 
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Inner  Radius 
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No.  of  Subrings 


Table  1 


Photodetector  diode  dimensions 


%  obscuration 

A 

B 

C 

D 

65.0 

1.1370 

0.1479 

0.9925 

0.01? 

79.4 

1.2597 

0.3005 

1.0250 

0.042 

87.8 

1.4425 

0.4855 

1.0575 

0.06^ 

92.5 

1.7094 

0.7433 

1.1010 

0.082 

97.5 

2.7160 

1.5244 

1.1850 

0.132 

Table  2  Correction  equation  coefficients  for  the 
Log-Normal  model. 


Liquid  volume  per  injection 
Injection  pump  speed 
Peak  injection  pressure 
Nozzle  orifice  diameter 


594.6  mm 
170  R.P.M. 
312  bars 
0.38  mm 
6.0  cm 


Distance  downstream  of  the  Nozzle 
Horizontal  distance  from  spray  axis 


0.8  cm 
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FIELD. 


PARTICLES 


HG'JRL  1.  OPTICAL  ARRANGED  NT  OF  FOURIER  TRANSFORM  LENS  TO  OBTAIN  THE  SIZE 
DISTRIBUTION  OF  SOLID  OR  LIQUID  PARTICLES. 
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FRACTION  OF  TOTAL  FNEROY  WITHIN  CIRCLE  OF  RADIUS 
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